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Output  Feedback  Tracking  Control 
of  an  Underactuated  Quad-Rotor  UAV* 

DongBin  Lee1,  Timothy  Burg1,  Bin  Xian2,  and  Darren  Dawson1 
September  25,  2006  (Updated) 


Abstract 

This  paper  proposes  a  new  controller  for  an  underactuated  quad-rotor  family  of  small-scale 
unmanned  aerial  vehicles  (UAVs)  using  output  feedback  (OFB).  Specifically,  an  observer  is 
designed  to  estimate  the  velocities  and  an  output  feedback  controller  is  designed  for  a  nonlinear 
UAV  system  in  which  only  position  and  angles  are  measurable.  The  design  is  performed  via  a 
Lyapunov  type  analysis.  A  semi-global  uniformly  ultimate  bounded  (SGUUB)  tracking  result 
is  achieved.  Simulation  results  are  shown  to  demonstrate  the  proposed  approach. 

Keywords:  Output  Feedback,  Observer,  Lyapunov,  Nonlinear,  Quad-rotor,  UAV,  Under¬ 
actuated 


1  Introduction 

The  potential  for  unmanned  aerial  vehicles  (UAVs)  in  applications  as  diverse  as  fire  fighting,  emer¬ 
gency  response,  military  and  civilian  surveillance,  crop  monitoring,  and  geographical  registration 
has  been  well  established.  Many  research  groups  have  provided  convincing  demonstrations  of  the 
utility  of  UAVs  in  these  applications.  However,  there  is  still  a  large  chasm  between  the  anticipated 
“tool  of  the  future”  and  currently  available  systems.  The  commercial  and  military  use  of  UAVs  is 
predicated  on  the  ability  of  such  vehicles  to  perform  new,  safer,  or  more  cost  effective  tasks  than 
traditional  manned  aircraft.  Until  recently,  this  has  been  more  of  a  question  than  a  statement; 
however,  recent  advances  in  aerial  vehicle  construction,  sensors,  digital  electronics,  control  design 
have  seen  a  rapid  increase  in  UAV  applications. 

Aerial  vehicle  construction  should  be  considered  as  an  important  factor  in  UAV  acceptance  and 
use.  Materials  such  as  carbon  fiber  can  be  used  to  reduce  weight  and  improve  robustness,  both 
critical  parameters  in  any  aerial  application.  Improved  manufacturing  techniques  are  capable  of 
producing  small,  complex,  precise  parts  at  a  reasonable  price  and  new  battery  technologies  have 
made  electric  hovering  craft  more  feasible.  One  of  the  interesting  small  aerial  vehicles  that  seems  to 
have  benefited  from  these  developments  is  the  quad- rotor  UAV  depicted  in  Figure  1.  The  quad-rotor 
consists  of  four  independently  driven  rotating  blades  that  can  provide  lift,  in  the  vertical  direction. 
The  vehicle  moves  in  other  directions  by  creating  a  mismatch  between  rotor  speeds,  and  hence,  this 
configuration  can  produce  torques  about  the  roll,  pitch,  and  yaw  axes.  The  basic  concept  for  the 
quad-rotor  dates  back  to  1907;  some  notes  on  the  history  of  the  quad-rotor  and  related  references 
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can  be  found  in  [1],  With  this  as  a  backdrop,  the  focus  of  the  work  presented  here  will  be  the  small 
quad-rotor  family  of  aerial  vehicles.  The  discussion  will  be  limited  to  vehicles  with  less  than  0.5kg 
payload.  This  weight  restriction  means  that  certain  technologies  that  may  make  sense  for  larger, 
more  expensive  aircraft  may  not  apply  to  this  class  of  aircraft. 

Position,  velocity,  and  acceleration  sensing  issues  will  begin  the  discussion  and  will  serve  as  par¬ 
tial  motivation  for  the  choice  of  output  feedback  control  design.  From  the  schematic  representation 
of  a  quad-rotor  UAV  shown  in  Figure  1,  it  can  be  seen  that  the  arbitrarily  positioned  aircraft  is 
fully  located  and  oriented  using  three  translational  positions  (x,  y,  and  z)  and  three  angles  (roll  - 
0,  pitch  -  0.  and  yaw  -  U).  Measurement  of  these  six  positions,  the  six  first  derivatives  (velocity  and 
angular  velocity),  and  the  six  second  derivatives  (acceleration  and  angular  acceleration)  has  proven 
challenging.  Measuring  angular  and  translational  quantities  each  present  different  challenges.  The 
angular  velocity  is  perhaps  the  most  accessible  angular  measurement.  Angular  velocity  can  be 
sensed  using  a  mechanical  gyroscope.  This  approach  has  historically  yielded  good  results  and  can 
be  scaled  down  in  a  cost  effective  manner.  The  major  drawback  is  the  moving  mechanical  parts 
may  add  additional  weight  and  reduce  reliability.  Piezoelectric  gyros  have  provided  an  alternative 
that  does  not  have  moving  parts.  More  reliable  gyroscopes  such  as  the  laser  ring  oscillator  are  not 
available  for  this  application.  Angular  position  may  be  more  difficult  to  measure  directly;  however, 
devices  such  as  magnetic  compass,  magnetometers,  tilt  sensors,  optical  horizon  sensors  may  provide 
estimates  of  position  for  the  roll  or  tilts  axes,  but  typically  have  low  bandwidth  and  poor  accuracy. 
There  are  several  technology  options  for  building  accelerometers  to  measure  the  angular  acceleration 
including  MEMs  and  piezoelectric  devices. 

Translational  position  can  be  directly  measured  with  global  positioning  system  (GPS)  based 
systems.  Enhancements,  such  as  DGPS,  are  required  to  achieve  improved  accuracy.  Direct  sensing 
of  the  linear  velocities  is  more  difficult;  an  anemometer  can  be  used  to  make  indirect  measurements 
that-  are  not  necessarily  relative  to  a  fixed  inertial  frame.  An  emerging  alternative  to  the  above 
sensors  is  to  use  vision  systems  to  measure  positions  or  velocities.  These  camera  systems  may  be 
ground-based  to  monitor  a  UAV  in  a  fixed  area  or  may  be  vehicle  based  and  used  to  estimate  changes 
in  scenery.  Finally,  it  is  often  difficult  to  convert  between  data  types  with  standard  mathematical 
operations.  For  example,  a  backwards  difference  estimate  of  velocity  from  position  information  may 
lead  to  a  noisy  velocity  signal,  and  integration  of  the  velocity  signal  to  obtain  position  information 
can  rapidly  accumulate  error.  If  a  trend  were  to  be  predicted  based  on  review  of  literature,  it  would 
be  that  angular  and  linear  positions  will  be  more  easily  attained  as  technology  evolves.  This  view  is 
demonstrated  in  [2]  where  only  a  single  GPS  sensor  is  used  to  measure  vehicle  position  and  velocity 
for  a  small  plane.  If  this  trend  is  true,  then  a  controller  that  uses  only  position  information  is  well 
motivated. 

The  generally  accepted  end-goal  that  a  vehicle  would  autonomously  take-off,  fly  to  a  mission 
site,  perform  a  mission,  return,  and  land  creates  a  daunting  challenge.  One  of  the  most  fundamental 
components  of  this  challenge  is  to  ensure  that  the  craft  can  move  to  or  hold  a  desired  position  and 
orientation.  Specifically,  as  shown  in  Figure  1,  the  aircraft  must  be  able  to  move  from  a  current 
location  to  a  new  desired  position  (denoted  by  the  triple  xd,  Vd-  zd)  and  achieve  a  new  orientation 
(denoted  by  the  angles  (j)d,  9d,  ripd  ).  This  low-level  control  objective,  as  it  often  called,  is  embedded 
at  the  center  of  high-level  control  objectives  such  as  path  planning,  target  tracking,  or  coordination 
with  other  crafts.  Design  of  the  low-level  control  represents  the  point  at  which  the  peculiarities  of 
the  mult-i-bladed  UAV  system  must  be  addressed;  specifically,  nonlinearities  of  the  system  dynamics 
and  the  fundamental  fact  that  the  system  is  under-actuated.  An  under-actuated  system  is  especially 
challenging  to  control  since  it  has  less  control  inputs  than  degrees  of  freedom,  i.e.  it  has  degrees  of 
freedom  that  cannot  be  directly  actuated.  In  order  to  achieve  high  overall  performance,  one  must 
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address  the  low-level  control  problem.  The  quad-rotor  UAV  has  six-degrees  of  freedom;  the  three 
translational  directions  and  the  three  rotational  angles;  however,  there  are  only  four  control  inputs; 
the  z-axis  thrust  and  the  three  rotational  torques.  The  quad-rotor  UAV  problem  is  sufficiently 
challenging  such  that  many  researchers  have  proposed  control  solutions  based  on  a  variety  techniques 
that,  accentuate  and  address  different  aspects  of  the  control  problem.  Solving  the  underactuated 
quad-rotor  problem  provides  a  choice  of  which  degrees-of-freedom  will  be  controlled.  For  example  in 
[1] ,  the  authors  use  feedback  linearization  to  explicitly  control  of  the  roll,  pitch,  and  yaw  angles  and 
the  height.  Translational  positions  are  then  implicitly  controlled  by  specification  of  the  trajectories 
for  the  controlled  axes.  More  typically,  the  translational  control  problem  is  directly  addressed  along 
with  yaw  angle  control.  For  example  in  [4],  the  authors  use  a  nested  saturation  algorithm.  In  [3],  a 
backstepping  approach  to  control  the  quad-rotor  based  on  a  model  of  the  specific  dynamics  of  the 
X4  flyer  is  used;  this  model  is  more  complicated  than  that  used  in  the  other  model-based  designs  as 
it  includes  additional  terms  such  as  gyroscopic  effects  of  the  rotating  blades.  In  [5],  the  authors  use 
a  quaternion-based  feedback  control  scheme  for  exponential  attitude  stabilization  of  a  quad-rotor 
aircraft.  The  work  given  in  [6],  [7],  and  [8]  are  representative  of  vision  based  applications. 

In  this  paper,  a  tracking  controller  is  designed  for  the  nonlinear  dynamic  model  of  the  quad-rotor 
helicopter  which  uses  only  output  feedback;  that  is,  the  controller  operates  using  only  position  and 
attitude  measurements.  To  appreciate  the  control  design  problem  it  is  useful  to  consider  the  quad- 
rotor  as  a  set  of  two  coupled  dynamic  subsystems.  The  first  subsystem  contains  the  translational 
dynamics  of  the  craft  and  has  a  single  input  along  the  z-axis  direction  and  the  second  subsystem  con¬ 
tains  the  rotational  dynamics  of  the  craft  and  contains  a  torque  input  to  each  of  the  three  rotational 
directions;  thus,  the  translational  subsystem  is  inherently  underactuated.  The  systems  are  coupled 
via  the  fact  that  the  rotational  velocities  appear  in  the  Coriolis  force  that  acts  on  the  translational 
subsystem  -  this  coupling  is  critical  as  it  provides  the  basis  for  the  backstepping  approach.  The 
choice  of  control  objectives  can  simplify  or  complicate  the  control  design  approach,  for  example, 
[Park]  sought  to  control  z-axis  position  and  the  three  angles  and  thus  the  control  inputs  are  already 
acting  at  the  point  of  interest  and  a  PID  control  was  directly  applied  (neglecting  nonlinearities). 
The  choice  here  to  control  the  three  linear  translations  and  the  yaw  angle  requires  that,  some  of  the 
torque  inputs  be  “redirected”  in  order  to  achieve  the  translational  tracking  objectives.  This  work 
builds  on  previous  backstepping  approaches  for  injecting  additional  control  action  into  the  transla¬ 
tional  dynamics  of  the  quad-rotor  but  is  greatly  complicated  by  the  co-design  of  a  velocity  observer. 
Velocity  estimation  is  a  well  known  problem  in  the  robotics  literature.  In  [9],  de  Queiroz  et.  al.  make 
a  systematic  presentation  of  the  observed  integrator  backstepping  technique  to  design  an  observer 
for  joint,  velocities  in  an  n-link  robot.  In  [9],  it  is  shown  that,  the  observer  and  controller  must, 
actually  be  designed  concurrently  via  Lyapunov  stability  arguments  in  order  to  ensure  a  stability 
result.  -  a  semi-global,  exponential  convergence  of  velocity  estimation  error  and  position  tracking 
is  shown  for  the  rigid  link  robot..  The  reader  is  referred  to  this  work  to  understand  the  observer 
and  controller  design  and  as  a  source  of  background  for  this  technique.  An  additional  reference  for 
the  observer  design  is  provided  to  [10]  where  this  same  observer  design  is  used  but  a  subtlety  of 
ensuring  that,  all  signals  in  the  implementable  form  of  the  observer  are  bounded.  The  payoff  from 
this  approach  is  the  mathematical  assertion  of  semi-global,  uniformly  ultimately  bounded  tracking 
while  compensating  for  system  nonlinearities,  estimation  error,  and  the  perturbations  that  result, 
from  the  backstepping  technique. 

The  paper  is  organized  as  follows.  In  Section  2,  a  well  known  model  of  the  quad-rotor  vehicle 
is  presented.  The  assumptions  and  properties  of  this  model  are  included.  In  Section  3,  a  velocity 
observer  for  output,  feedback  tracking  control  is  designed  to  estimate  the  velocities  and  the  OFB 
controller  using  this  velocity  estimate  is  derived.  Stability  analyses  on  observer,  controller,  and 
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composite  system  are  considered  mathematically  in  Section  4  followed  by  Theorem  1  and  Remark 
1  which  show  the  analysis  result  and  the  boundededness  of  all  signals.  Simulation  results  demon¬ 
strating  the  performance  are  presented  in  Section  5.  Detailed  developments  and  stability  proofs  are 
deferred  to  Appendices  followed  by  a  conclusion. 


2  System  Modeling 

2.1  Quad-Rotor  Aerial  Vehicle  Model 

The  quad- rotor  unmanned  aerial  vehicle  is  shown  in  Figure  1  where  it  can  be  seen  that  the  body- 
fixed  reference  frame,  B,  moves  relative  to  a  fixed  inertial  frame,  I.  The  translational  and  rotational 
dynamic  equations  of  motion  in  the  body-fixed  reference  frame  are  [6] 

ml3  03x3 
03x3  J 


where  v(t)  G  5ft3  denotes  the  linear  velocity,  uj(t)  =  [ojx,ujy1ujz]T  G  5ft3  represents  the  angular 
velocity,  m  G’ft1  is  the  known  mass  of  the  quad-rotor,  J  G  5ft3x3  denotes  a  positive  definite  inertia 
matrix,  G(R)  G  ft3  is  a,  gravity  vector  described  below,  and  Ni(v),  N2(u)  G  !ft3  are  the  nonlinear 
aerodynamic  damping  interactions.  The  input  u\[t)  G  ’ft1  provides  lifting  force  in  the  z-direction 
and  U‘i{t)  G  ’ft3  creates  rotation  torque  in  the  roll,  pitch,  and  yaw  directions.  The  specific  form 
of  the  quad-rotor  links  the  inputs  to  the  dynamics  via  Bi  =  [0,0, 1]T  G  ft3  and  Ii2  =  I3  G  5ft3x3. 
Additionally,  I3  is  a  3x3  identity  matrix,  03xi  G  ’ft3  represents  a  3x1  zero  vector  and  03x3  G  ’ft3x3 
represents  a  3x3  zero  matrix,  and  S(-)  G  ’ft3x3  is  a  general  form  of  skew- symmetric  matrix  as  follows 
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It  will  be  necessary  to  relate  the  inertial  frame  coordinates  to  the  body-fixed  frame  coordinates  as 
seen  in  Figure  1.  The  quaternion  based  rotation  matrix  R(q)  G  5ft3x3  that  translates  a  body-fixed 
frame  vector  into  inertial  coordinates  is  calculated  from  the  following  form 


R(q)  =  ( qI  ~  QvVv)! 3  +  2 qvql  -  2 qQS(qv) 


(3) 


where  q(t)  =  [q0,qv]T  ^  represents  the  unit  quaternion,  qQ(t )  G  5ft1,  qv(t )  G  5ft3  are  subject 
to  the  constraint  q%  +  qj,qv  =  1  [11].  The  unit  quaternion  can  be  generated  from  u(t)  by  the 
relationship  [12] 
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The  angular  velocity  transformation  matrix  based  on  the  Euler  angles  [12],  denoted  as 
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Figure  1:  Quad- Rotor  UAV  Coordinates 

is  used  to  relate  the  rate  of  change  to  the  angular  velocities  in  the  body-fixed  frame  by  the  following 
expression 

0  =  T(Q)u  (6) 

where  the  Euler  angles,  0(f)  =  [</>,  9,'ip]T  G  Ji3,  describe  the  orientation  of  the  body-fixed  frame 
relative  to  the  inertial  frame  and  in  which  c-  =  cos  (■ ),  s ■  =  sin  ( • ) ,  t-  =  tan  (■)  are  used.  Note 
that  we  use  the  quaternion  based  translational  rotation  matrix,  R(q),  in  order  to  avoid  the  matrix 
singularity  associated  with  the  Euler  angles  representation  where  possible;  however,  the  Euler  angles 
based  angular  transformation  matrix  must  be  used  to  facilitate  access  to  the  yaw  angle  for  tracking 
control.  A  gravity  vector  is  denoted  as 

G(R)  =  mgRT(q)Ez  G  5i3  (7) 

where  g  G  Ji1  denotes  the  gravitational  acceleration  due  to  the  gravity  and  Ez  =  [0,0,  l]T  denotes 
the  unit  vector  in  the  coordinates  of  the  inertial  frame.  G(R)  is  represented  in  the  body- fixed  frame 
as  the  multiplication  by  RT (q)  would  indicate. 

The  translational  and  rotational  kinematic  equations  in  the  body-fixed  reference  frame  are  given 

by 
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_wj  L  °3x3  T-\Q) 

where  p(t)  e  T3  contains  the  position  of  the  body-fixed  reference  frame  relative  to  the  inertial 
frame,  and  its  derivative  pit)  G  T,!  represents  the  translational  velocity  in  the  inertial  frame.  The 
transformation  in  (8)  can  be  reduced  to  the  simplified  form 

$(v,u)=D(R,T)x  (9) 

where 
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The  dynamics  in  (1)  can  be  compacted  and  transformed  into  the  inertial  frame  as  shown  in  Appendix 
A  to  yield  the  dynamic  model 

M*(R,  T)x  =  N*(R,  T,  x)x  +  h*(R ,  T,  x )  +  G*(R )  +  B*(R,  T)U  (12) 

where  M*(R:T )  E  '!RGx6  denotes  the  inertia  matrix,  N*(R,  T.  x)  E  sRGxG  is  a  Centrifugal/ Coriolis 
force  matrix,  h*(R,  T,  x)  G  T6  is  a  hydrodynamic  damping  term,  G*(R)  G  'h6  is  a  gravity  term,  and 
B*(R,T)  G  Ji6x4  represents  the  input  matrix.  All  are  explicitly  defined  as  (143)  in  Appendix  A. 

2.2  Model  Properties  and  Assumptions 

2.2.1  Model  Properties 

The  dynamic  system  given  in  (12)  satisfies  the  following  properties. 

PI:  The  inertia  matrix  M*(R,T)  and  Centrifugal/ Coriolis  matrix  N*(R,T,x)  in  (12)  satisfy  the 
following  skew-symmetric  property  [9] 

£t  (J-t  (. M*{R ,  T))  +  2 N*(R,  T,  x)^  =  0,  V  £  G  d£6, 

which  is  proven  in  Appendix  B. 

P2:  The  inertia  matrix  M*(R,  T)  can  be  upper  and  lower  bounded  in  the  following  form 

mr  lief  <  tTM*(R,  T)£  <  m2  ||£||2  ,  V  £  G  K6 

where  m i,  m2  G  k1  are  positive  constants. 

The  rotation  matrix  R(q)  of  (3)  where  R(q )  =  R(Q)  and  the  skew- symmetric  matrix  S(u)  of 
(2)  satisfy  the  following  properties  [12]. 

P3:  R  =  RS{u),  Rt  =  -S(co)RT 

P4:  Rt  =  R~\  RtR  =  RRt  =  I3 

P5:  ST(()  =  -S(0,  s(()6  =  -sm,  V{,  6  e  K3 

2.2.2  Model  Assumptions 

The  following  assumptions  are  made  regarding  specific  components  of  the  dynamic  model.  Note 
that  this  assumptions  can  be  shown  true  for  specific  models  but  not  in  general. 

7f 

Al:  With  regard  to  angular  transformation  matrix,  T(0)  introduced  in  (5),  we  assume  that  0  ^  ±  — 

2 

or  T-1(0)  exists  and  ||T(0)||i  is  bounded  (i.e.,  ||T,(0)||i  <  e\  where  E\  G  5R1  is  a  positive 

constant). 
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A2:  The  nonlinear  damping  term  Ni(v)  G  P3,  N2(u>)  G  P3  introduced  in  (1)  can  be  replaced  by  the 
linearly  parameterized  form  Y1(v)9i  =  Ni(v),  Y2(lu)92  =  N2(u)  where  Y1(v)  G  5£3xZ,  Y2(cu)  G 
3?3xm  are  known  regression  matrices  and  9\  G  iR1,  92  G  Rm  are  known  parameter  vectors.  The 
dimension  of  l  and  m  are  related  to  the  specific  model  for  Ni(v),  N2(uj).  Additionally,  it  is 
assumed  that  the  first  nonlinear  parameterized  term  can  be  upper  bounded  as 

HUMMI  <fc4IMI,  where  £c4  G  S'  is  a  positive  constnat. 

A3:  The  nonlinear  term  h*(-)  G  T6  in  (12)  can  be  linearly  parameterized  in  terms  of  x(t)  as  follows 

h*(R,T,x)  =  Y*(R,T)x  (13) 

in  which  Y*(R.  T )  G  TGx?)  and  Y*(R.  T )  can  be  upper  bounded  in  the  following  form 

||F*(.R, T)||  <  £c0,  where  £c0  G  3?1  is  a  positive  constant.  (14) 

3  Output  Feedback  Tracking  Control 

The  goal  of  the  tracking  controller  is  to  force  the  aerial  vehicle  track  a  desired  trajectory.  As 
discussed  in  the  modeling  section,  the  quad-rotor  aerial  vehicle  is  under-actuated,  and  hence,  a 
decision  must  be  made  as  to  which  degrees  of  freedom  are  to  be  controlled.  First  the  choice  was 
made  to  control  the  translational  position,  p(t)  G  P3  in  the  inertial  reference  frame,  along  with 
yaw,  Y{t)  G  P1  in  the  inertial  reference  frame.  The  translational  position  p(t)  and  the  angular 
position  (~)(t)  G  P3  are  the  only  measurable  states,  other  states  such  as  the  translational  velocity 
v(t)  G  P3  and  the  angular  velocity  oj(t)  G  P3  are  not  measurable  and  cannot  be  included  in  the 
controller  design.  Here  we  assume  that  the  desired  trajectories  and  their  up  to  third  derivatives  are 
all  bounded;  i.e.,  pd(t),  pd(t ),  pd(t),  and  Pd  (t)  G  Too  and  ipd(t),i>d(t ),  and  ipd(t)  G 

3.0.3  Full  State  Feedback  (FSFB)  Error  Systems  to  Motivate  the  Structure  of  the 
Output  Feedback  Control 

In  order  to  demonstrate  the  approach  to  designing  an  OFB  controller  the  same  approach  is  demon¬ 
strated  for  the  less  complicated  FSFB.  The  position  tracking  error,  denoted  as  ep(t)  G  T3,  expressed 
in  the  body-fixed  frame  is  defined  as  the  transformed  difference  between  the  inertial-frame  based 
position  and  the  inertial-frame  based  desired  position,  denoted  as  pd(t)  G  T3,  in  the  manner 

ep  =  RT(p-pd).  (15) 

The  position  tracking  error  rate,  ep  (t)  G  T3,  is  obtained  by  taking  the  time  derivative  of  (15), 

ep  =  RT(p-  pd )  +  RTp  -  RT pd)  ( 16) 

and  after  substituting  for  R1  =  —S{uj)Rt  from  P3,  and  pit)  =  Rv  from  (8),  and  using  RTR  =  J3 
in  P4  we  have 

ep  =  -S(uj)Rt(p  -pd)  +  v  -  RTpd-  (17) 

We  now  use  the  definition  of  ep(t)  in  (15)  to  collect  terms  in  (17)  and  then  the  term  ^ RT pd(t )  is 
added  and  subtracted  to  facilitate  the  introduction  of  the  term  ev(t)  G  P3  as  follows 

ep  =  -S(u)ep  +  —  (mv  -  RTpd )  +  —  RTpd  -  RTpd.  (18) 
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From  the  equation  in  (18)  the  translational  velocity  tracking  error,  ev(t)  G  ft3,  in  the  body-fixed 
frame  is  defined  as 

ev  =  mv  —  Rt pd  (19) 

where  pd(t)  G  5ft3  is  the  desired  translational  velocity.  The  final  form  of  the  open- loop  position 
tracking  error  dynamics  in  full  state  feedback  system  is  obtained  from  (18),  (19)  as  follows 

ep  =  —S(uj)ep  +  —ev  +  (— —  l)RTpd.  (20) 

mm 

After  taking  the  time  derivative  of  ev(t)  in  (19),  then  substituting  for  mv(t)  from  (1)  and  for 
RT (•)  =  —S(uj)Rt  from  P3,  we  have 

ev  =  -mS(u)v  +  Y1(v)9i  +  G(R)  +  BiUx  +  S(cu)RTpd  -  RT  Pd  (21) 

where  A2  was  used  to  replace  Ni(v).  After  collecting  the  terms  in  (21)  and  applying  the  ev(t) 
definition  in  (19),  we  have  the  velocity  error  rate  as 

ev  =  -S(u)ev  +  G(R)  +  Yl(v)0l  -  RT  Pd  +BlUl.  (22) 

The  yaw  angle  tracking  error,  e$(t)  G  5ft1,  is  defined  in  the  inertial  coordinate  system  as 

e^  =  ip-  i)d.  (23) 

The  goal  in  the  control  development  will  be  to  ensure  that  e^(t)  and  ep(t)  are  driven  to  small  values; 
that-  is,  to  ensure  the  control  objectives  are  met.  The  yaw  angle  rate  error  system  is  derived  by 
taking  the  time  derivative  of  (23)  as  follows 

e^  =  'i/j-'ijjd  =  Tz(Q)u-uzd  (24) 

where  the  term  TZ(Q)  G  5ft 1x3  is  the  third  row  of  T(0)  from  (5).  Note  that  Tz(Q)lu  =  ^  in  (6)  and 
rjjzd  =  '0d  where  uzd(t )  is  the  desired  yaw  angle  rate  in  the  body-fixed  frame.  In  order  to  further 
develop  the  control  design,  the  filtered  position  tracking  error  signal  rp(t)  G  5ft3  is  defined  in  the 
following  manner  [6] 

rp  =  ev  +  aep  +  8  (25) 

where  a  G  5ft1  is  a  positive  constant  and  8  =  [  0  0  <53  ]T  G  5ft3  is  a  constant  design  vector  in  which 
hs  G  5ft1  is  a  scalar  constant.  The  filtered  position  tracking  error  can  be  combined  with  the  yaw 
tracking  error  to  create  a  composite  tracking  error  r(t)  G  5ft4  in  the  manner 


r  = 


if; 


(26) 


The  filtered  tracking  error  dynamics  can  be  found  by  first  differentiating  (26)  to  yield 


(V> - 

T 

1  V 

CV  -\~OL  Gp  1 

/  - 

i 

Cb- 

•0) 

_ 1 

(27) 


The  filtered  position  tracking  error  rate,  fp(t ),  are  obtained  by  substituting  (20)  and  (22)  to  yield 


Q,  1 

fp  =  —S{iS){ev  +  otep  +  8  —  8)  -\ - ev  +  «( - 1  )RTpd  —  RT  Pd  +G(i?)  +  ki(u)di  +  B\U\  (28) 

mm 
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where  the  term  S(u)6  has  been  added  and  subtracted  to  facilitate  introduction  of  rp(t)  G  ’ft3  on  the 
right-hand  side  as  shown  below 

A/  1 

rp  =  —S{uj^rp  +  - — ev  +  a(- - l)RT,pd  —  RT  Pd  -\-G(R)  +  Y\(v)9i  +  +  B\U\\ .  (29) 

It  is  now  a  straightforward  matter  to  substitute  from  (24)  and  (29)  into  (27)  to  yield  the  open- loop 
filtered  tracking  error  dynamics  in  the  following  form 


— S(u)rp  +  ^ ev  +  a(^  —  1  )RTpd  —  RT  Pd  +G(R)  +  Yi{v)6\ 

_|_ 

—S(8)u  +  7?i»i 

—uzd 

1 

Tz(Q)uj 

where  S(£)8  =  —S(8)£  in  P5  was  used  to  modify  the  S(u)6  term.  The  last  square  bracketed  term  in 
(30)  highlights  the  location  where  the  control  input  will  be  eventually  designed.  The  control  input, 
Ui(t)  can  be  designed  by  introducing  the  auxiliary  signal  z(t)  G  T3  in  order  to  inject  an  auxiliary 
control  signal  u\(t)  into  the  translational  dynamics  from  the  rotational  dynamics,  uj(t)  as 


z  =  u  — 

where  Bz  =  [I3,  03xi]  G  '!R3x4  and  the  control  input 

u\  =  [  0  0 


Bzui 

(31) 

is  defined  by 

0  1]  iii 

(32) 

Then  we  have 


r  = 


-S(uj)rp  +  ^ev  +  a(±  -  l)RTpd 


RTPd  +G(R)  +  Y]_(v)9i  +  ^(e 
-uzd 


p 


+Bpui+Bp 


(33) 


where  the  procedure  is  detailed  in  the  Appendix  C  titled  “Development  and  Derivative  of  Con¬ 
trol  Signal  the  term  ^ep(t)  was  added  and  subtracted  in  (30)  to  facilitate  further  control 

development,  and  Bp(-)  G  '!R4x4  is  defined  as 


B 


-S{6)  Bx 
TZ(Q)  0 


(34) 


3.1  Observer  Design 


The  next  step  in  the  control  development  is  to  address  the  problem  that  the  velocity  of  the  aerial 
vehicle,  x(t)  is  not  directly  measurable.  In  the  inertial  coordinate  system  the  vehicle  velocity  and 
vehicle  angular  velocity,  x(t).  cannot  be  obtained  from  sensor  readings.  The  well  known  method 
of  circumventing  this  problem  is  to  create  an  estimate  of  the  unmeasurable  states  for  use  in  the 
feedback  control  law.  The  estimated  linear  and  angular  velocities  v(t),  u(t)  G  T3  are  introduced 
using  (9)  as 

RT  (q)  03x3 

o3x 3  t~\q) 


V 

Cj 

P 

J 

.  © 

(35) 


T  T 


The  equation  (35)  can  be  rewritten  by  introducing  $(D,tD)  =  [vT ,u> 
to  yield 

4,(h,d))  =  D(R,  T)  x  . 


■  T 


•  T 


and  x  (t)  =  \p  ,0  }7 
(36) 
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The  goal  of  estimating  velocity  is  achieved  by  using  an  observer  similar  to  [9]  as  follows 


x=  y  +  koix  (37) 

where  the  initial  condition  y( 0)  =  — Tcoi^(O) ,  koi  G  3?1  is  a  positive  gain  and  x(t)  =  [ pT ,  0T]  G  J£6  is 
defined  as 

x  =  x  —  x  (38) 

where  x(t)  =  [ pT ,  0T]T  G  3?6  denotes  a  position  estimate  vector.  The  auxiliary  signal  y(t)  G  J£6 
introduced  in  (37)  is  updated  according  to 

y  =  M*(R,T)~ 1  (N*(R,T,x0)x0  +  h*{R,T,x0)  +  G*(R)  +  B*{R,T)U  +  k02x)  +  k03x  (39) 

where  k0 2,  k03  G  9?1  are  positive  constants  and  an  auxiliary  signal  x0(t)  G  5£6  is  defined  as 

x0  =x  —/3x  (40) 

where  (3  G  3?1  is  a  positive  constant.  The  equation  (39)  was  obtained  by  substituting  x0(t)  for  x(t) 
into  the  dynamic  model  in  (12)  and  adding  the  x-terms  to  promote  stability  of  the  observer.  The 
basic  form  of  the  internal  observer  signal  y(t)  G  J£6  comes  from  substituting  x0(t)  into  the  modeling 
equation  of  (12);  that  is,  the  knowledge  of  the  system  dynamics  and  parameters  are  exploited  to 
produce  the  velocity  estimate.  In  order  to  simplify  the  Lyapunov  analysis,  the  following  filtered 
observer  error  signal  is  defined  simlar  to  [10]  as 

s  =  x  —  xa  =x  +/3x  G  5i6  (41) 


where  (40)  was  used.  This  estimate  must  be  faithful  to  the  actual  velocity  signal  in  the  sense  that 
the  velocity  estimate,  x  (■ t ),  will  produce  a  small  value  for  the  velocity  estimation  error,  x  (t)  G  5i6, 
defined  as 

x=x  —  x  (42) 

where  the  estimated  velocity  errors  can  be  represented  by  defining 


4>(fi,a;) 


=  D(R,T )  x=  D(R,T ) 


V 

0 


(43) 


and 


p=  p —  p  and  0= 


0-  0  . 


(44) 


The  motivation  for  the  form  of  the  observer  and  the  proof  that  the  observer  will  produce  a  proper 
estimate  in  the  observation  error  dynamics  is  developed  below.  Taking  the  time  derivative  of  (37) 
yields 

x=  y  +  koi  x  .  (45) 

Multiplying  (45)  by  M*(R,T):  and  then  substituting  from  (39)  for  M*(R,T)y(t)  yields 


M*(R,T)  x 


N*(R,  T,  x0)x0  +  h*(R,T,x0)  +  G*(R)  +  B*(R,T)U 
+k02x  +  k03M*(R ,  T)x  +  k01M*(R ,  T)  x  . 


(46) 
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After  subtracting  (46)  from  (12),  the  dynamics  of  x(t)  is  as  follows 

M*(R,T)'x  =  N*  (R,  T,  x)x  —  N*  (R,  T,  x0)x0  +  h*  (R,  T,  x)  —  h*  (R,  T,  xa) 

-k02x  -  k03M*(R,  T)x  -  k01M*(R,  T)  x  .  (47) 

On  the  other  hand,  taking  time  derivative  of  (41),  and  multiplying  M*(R,  T )  yields 

M*(R,T)s  =  M*(R,T )  x  T)  x  .  (48) 

Substituting  (47)  for  M*(R,T )  x  ( t )  and  arranging  the  terms  yields 

M*(R,T)s  =  N*(R,T,x)x-N*(R,T,x0)x0  +  h*(R,T,x)  -h*(R,T,  x0) 

-[fc02  +  k03M*(R,  T)]x  +  [P  -  km\M*(R,  T )  £  .  (49) 

Introducing  the  new  positive  constants 

*03  =  *03/3  and  k01  =  (3  +  *03  (50) 

allows  the  last  two  terms  in  the  right  side  of  the  (49)  to  be  rewritten  as 

M*(R,T)s  =  N*(R,  T,  x)x  —  N*(R,  T,  x0)x0  +  h*(R,  T,  x)  —  h*(R,  T,  xQ) 

-k02x-k03M*(R,T)s  (51) 

where  (41)  was  used.  The  definition  of  s(t)  in  (41)  is  substituted  as  x(t)  =  s(t)  +  xD(t)  into  the 
outside  term  of  the  first  equation  to  the  right  side  in  (51)  and  then,  we  have 

M*(R,  T)s  =  N*(R,  T,  x)s  +  N*(R,  T,  x)xa  -  N*(R,  T,  x0)x0  +  h*(R,  T,  x)  -  h*(R ,  T,  xQ) 

-k02x  -  k03M*(R ,  T)s.  (52) 

Substituting  A3  yields 

M*(R,  T)s  =  N*(R,T,x)s+N*(R,T,x)xo-N*(yR,T,xo)xo+Y*(R1T)s-k02x-k03M*(yR,T)s.  (53) 

The  Centrifugal/Coriolis  terms,  N*  (R,  T,  x)x0—N*  (R,  T,  x0)x0  are  now  utilized  to  write  N*  (R,  T,  x0)s 
developed  in  the  Appendix  F.  Therefore,  the  velocity  observer  error  dynamics  yields 

M*(R,  T)s  =  N*(R,  T,  x)s  +  N*(R,  T,  xQ)s  +  Y*(R ,  T)s  -  k02x  -  k03M*(R,  T)s  (54) 

3.2  Output  Feedback  Control  Formulation 

Figure  2  shows  the  outline  of  the  output  feedback  tracking  system.  We  designed  the  observer  to 
estimate  the  velocities  v(t),  u(t )  which  are  not  measurable,  producing  x  (t)  and  then  closed- loop 
controller  using  backstepping  approach  based  on  the  Lyapunov  stability  which  will  enable  to  track 
the  desired  position  and  yaw  angle  is  designed  using  the  output  feedback  and  estimated  states  based 
on  the  error  dynamics  of  the  underactuated  UAV  system. 
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x=y+k^x,  y=f(RTJ,x0,U)+k& 


Figure  2:  Block  Diagram  of  OFB  Position  and  Yaw  Tracking  System 


3.2.1  Error  Definition  in  Output  Feedback  Control 

In  the  OFB  only  the  position  tracking  error  in  (15)  and  yaw  angle  tracking  error  in  (19)  are  same 
those  of  full  state  feedback  system.  The  definition  of  velocity  tracking  error  in  (19)  is  now  used  as  a 
guide  to  account  for  the  fact  that  the  velocity  is  not  measurable.  Specifically,  an  observed  velocity 
tracking  error  signal,  ev(t)  G  P3,  can  be  defined  by  substituting  the  estimated  linear  velocity,  v(t). 
into  (19)  as 

ev  =  mv  —  RTpd ■  (55) 

The  velocity  error  rate  equation  in  (22)  is  redefined  by  substituting  v(t)  and  ih(t)  to  yield  the 
observed  vlocity  rate  error  as 

ev=  —S(ih)ev  +  G(R)  +  Yi(v)9i  —  RT  Pd  -\-B\U\.  (56) 


In  a  similar  fashion,  ep  (t)  can  be  formed  from  (20)  by  substituting  Cj  and  ev 

ep=  ~S(uj)ep  +  —  ev  +  ( - 1  )RTpd-  (57) 

mm 

The  new  signal  f(t)  G  F4  has  been  introduced  to  represent  the  observed  filtered  tracking  error  and 
is  defined  by 


- 1 

e„  +  otep  +  6 

The  component  definition  of  r(t)  in  (58)  can  be  differentiated  to  yield 


rp 

ev  +cx  ep 

- 1 

■-9- 

1 

■-9- 

_ 1 

(58) 


(59) 


The  observed  filtered  position  error  rate,  fp(t)  is  found  by  substituting  (56)  and  (57)  into  the  top 
row  of  (59)  to  yield 


a  „  .1 

— ev  +  co¬ 
rn  TO 


1  )RTpd  +  G{R )  +  Yi(h)6h  —  Rt  Pd  +  —S(6)u  +  B\Ui] . 


(60) 
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The  rate  of  yaw  error  is  found  by  substituting  u(t)  into  (24) 


-ujzd  =  Tz(Q)uj  -  ujzd  (61) 

where  k'd(t)  =  uzd.  The  observed  tracking  error  rate  is  obtained  by  substituting  (60)  and  (61)  as 


r= 


S(u)rp  +  ^ ev  +  a(A  —  1  )RTpd  +  G(R)  +  Yi(v)6i  —  RT  Pd 

_|_ 

—S{8)u  +  B\Ui 

~UJzd 

1 

Tz(Q)u 

(62) 


The  equation  (62)  can  be  derived  by  following  the  procedure  which  is  shown  the  equations  from 
(30)  to  (34)  as 


r= 


—S(uj)rp  +  ^ ev  +  —  1  )RT,pd  +  G(R )  +  Y\{v)9\  —  RT  Vd 


-uzd 


Bpu  i  +  Bp 


0 


(63) 


where  the  auxiliary  signal  z(t)  G  $i3  is  introduced  in  the  same  way  in  (31)  by  substituting  u(t)  to 
denote 

z  =  lu  —  Bzui ,  (64) 

and  hence,  the  auxiliary  estimation  signal  z(t)  G  T3  can  be  defined  using  the  (31)  and  (64)  as 


Z  =  IjU  —  U)  =  UJ . 


(65) 


3.3  Controller  Formulation 

3.3.1  Translational  Input  as  a  Lifting  force 

The  filtered  tracking  error  dynamics  in  (33)  now  presents  the  opportunity  to  design  the  auxiliary 
control  input  v,\  (t).  It  can  be  seen  from  (33)  that  certain  terms,  those  containing  measurable  signals, 
can  be  directly  canceled  by  design  of  Ui(t).  Other  unmeasurable  terms,  e.g.,  ev(t),  rp(f),  Ti(u)6h, 
can  be  canceled  to  the  best  degree  possible  using  an  estimate  of  velocity.  Additionally,  a  term  of 
the  form  — r(t )  G  k4  will  be  required  to  promote  the  convergence  of  r(t)  to  zero.  With  these  three 
objectives  in  mind,  the  control  U\{t)  in  (33)  can  be  designed  based  on  equation  (63)  as 

~m^v  ~  a(m  ~  l)RTPd  -  G(R)  -  Y1{V)61  +  RTpd  -  ±ep 

uzd 


where  kr  =diag(/cri,  kri,  kri,  kr 2)  G  3?4x4  is  a  positive  constant  matrix.  For  stability  analysis  of  this 
control,  the  observed  filtered  tracking  error,  denoted  by  f(t)  G  T4,  can  be  described  by  subtracting 
(58)  from  (26)  as 


'  p 
0 

where  the  observed  filtered  position  tracking  estimation  error,  rp(t)  G  F3 
difference  between  rp(t )  in  (25)  and  rp(t)  in  (58)  and  can  be  shown  to  be 

rp  =  rp  —  rp  =  mv. 

The  difference  between  the  actual  velocity  error,  ev(t),  and  the  estimated 
defined  as 


(67) 

is  now  defined  as  the 

(68) 

velocity  error,  ev(t)  is 


e-v  =  ev  —  ev  =  m(y  —  v)  =  mv  =  rp 


(69) 
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where  v(t)  =  v  —  v  is  the  velocity  tracking  error.  The  equation  (66)  can  be  represented  by  substi¬ 
tuting  —r(t)  +  r(t)  for  —r(t)  and  substituting  —ev(t)  +  ev(t)  for  —ev(t)  in  the  following  form 


~mev  +  ~  a(m  ~  l)RTPd  ~  G{R)  ~  Y1(v)9i  +  RTpd  ~  ±ep  1 

Uzd  J 

(70) 

Finally,  the  closed-loop  filtered  tracking  error  dynamics  for  r(t)  is  formed  by  substituting  (70)  into 
(33)  to  yield 


Ui  =  -£>m  1  <  —krr  +  krr  + 


r=  —krr 


SHrp-^ 

0 


—ev  +  Yl91 

m  u  x  x 

o 


where  (67),  (69),  and  (31)  were  used,  and  Yi(-)  G  3?3x/  is  introduced  as  follows 

Y1(y)  =  Y1(y)-Yl(v). 


(72) 


3.3.2  Development  of  the  Torque  Input  using  a  Backstepping  Approach 

Examining  the  meaning  of  the  term  z(t)  introduced  in  (33)  should  help  crystallize  the  exposition  of 
the  control  design  approach  and  motivate  the  next  step.  The  definition  of  z(t )  in  (31)  quantifies  the 
closeness  of  the  control  term  U\(t)  to  the  angular  velocities  cn(f),  if  z(t)  is  zero  then  uj(t)  =  U\(t). 
The  implication  is  that  the  effect  of  rewriting  the  input  term  in  (33)  was  to  inject  the  signal  Ui{t) 
as  a  desired  input  to  the  rotational  dynamics  (a  backstepping  approach).  The  design  now  proceeds 
to  ensure  that  the  auxiliary  signal  z(t)  in  (31)  is  driven  to  a  small  value.  Taking  the  time  derivative 
of  z(t )  in  (31)  and  multiplying  by  the  inertia  matrix,  J,  yields 

Jz  =  JCj  —  JBZ  ui  .  (73) 

Substituting  the  second  equation  of  (1)  for  Jio(t)  into  (73)  produces 

Jz  =  S(Jlu)lu  +  N2(cl>)  +  B2u2  —  JBZ  ui  (74) 

where  it  can  be  viewed  as  backstepping  -v,\  (t)  through  the  integrator.  It  is  now  useful  to  group 
terms  in  equation  (74)  and  invoke  A2  for  the  parameterization  of  N2(oj)  as 

Jz  =  (^S(Juj)uj  +  Y2(u)92  —  JBZ  u^j  +  B2u2.  (75) 

The  following  assumption  is  made  for  (75): 

A4:  A  linear  parameterization  has  been  assumed 

Y3(p,  R,  v,  cu)93  =  S(Ju)uj  +  Y2(uj)92  -  JBZ  u\  (76) 

where  Y3(p,  R,v,cu)  G  T3x”  and  93  G  ft”  (i.e.,  there  are  n  known  parameters  in  93  and  n  is 
determined  by  the  specific  models  of  (76),  especially  Ni(v),  N2(lu)  and  (79)). 

Using  the  A4,  (75)  is  rewritten  as 

Jz  =  Y3(p,  R,v,uj)93  +  B2u2.  (77) 
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It  should  also  be  clear  that  ideally  the  control  input  u2(t)  would  be  designed  to  stabilize  the  z(t)~ 
dynamics  and  cancel  Y3(p,  R,v,u)9 3,  of  course  this  can  not  be  achieved  directly  because  both  of 
these  objectives  would  require  knowledge  of  unmeasurable  quantities.  The  estimated  velocities  will 
provide  the  best  opportunity  to  achieve  these  goals;  that  is,  in  a  similar  way  in  (77),  the  estimated 
velocities  v(t)  and  u(t)  are  substituted  into  (76)  to  create  the  estimate  Y3(p,  R,v,lu)  G  di3xn  given 
by 

Y3(p,  R,  V,  u)03  =  S(JCj)<2)  +  Y2(u)92  -  JBZ  u3  (78) 

where  Y2(i2j),  Y3(p,  R,v,u)  are  the  estimated  regression  matrices.  The  time  derivative  of  Ui(t)  can 
be  calculated  using  the  error  definition  in  (66)  as  follows 

(79) 

where  U(t)  is  from  the  parenthetical  terms  on  the  right  equation  (66)  and  the  time  derivative  of 
U  ( t )  is  defined  as  follows 


tU  = 

dt 


—kr  r  + 


■—  (Yi(,v)9d 

m  dt 


1  ■ 


m 


0 


S(u)  (  * RTpd  -  aRTpd  -  RT  Pd  +mgRTE 2 

Uzd, 


a  (1  -  i)  RT  Pd  +RT  Pi 


(80) 


where  (56),  (57),  (63),  and  G(R)  =  —S(cu)G(R)  are  utilized  and  the  time  derivative  of  Yi(v)$i  are 
explicitly  calculated  in  Appendix  E. 

The  control  input  u2(t)  G  di3  is  now  formulated  from  (77),  making  use  of  (78),  in  the  following 
form 

u2  =  B i1  ( -kzz  -  Y3(p,  R,  v,  u)63  -  Bj/r )  (81) 

where  the  first  term  is  a  linear  feedback  control,  the  last  term  is  added  to  cancel  a  crossing  term 
during  the  Lyapunov  stability  analysis,  and  Zi/;  ( • )  G  !R4x3  is  formed  from  the  first  three  columns  of 
B^-)  in  (34)  and  can  be  transposed  as  follows 

Bl  =  [-S(«)T.  77(6)T].  (82) 

After  substituting  (81)  into  (77),  —z(t)  +  z(t)  for  z(t),  and  using  u(t)  =  z(t )  in  (65),  we  have 

Jz  =  -kzz  +  kzG  +  Y393  -  B^r  +  B^f  (83) 


where  (67)  was  used  to  create  the  last  two  terms,  and  the  regression  estimation  error,  Y3(-)  G  T3xn, 
is  defined  as 

=  Y3(p,  R,  v,  u)  -  Y3(p,  R,  v,u).  (84) 

ST(0  =  in  P5  can  be  invoked  to  rewrite  the  matrix  -Bj(-)  in  (82)  as 

Bl  =  [S(6),  Tz(Q)t],  (85) 


and  hence,  we  have 


Bl f  =  |S(«),  rs(9)T] 


mv 

0 


=  mS(6)v. 


(86) 


After  substituting  (86)  into  (83),  we  have  the  final  form  for  the  closed- loop  system  as  shown  below 


Jz  =  —kzz  +  kzoj  +  Y393  —  BTr  +  mS(8)v. 


(87) 
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4  Stability  Analysis 

The  combination  of  velocity  observer  error  and  closed- loop  error  systems  given  by  (54),  (71),  and 
(87)  yields  the  following  stability  result  for  the  velocity  observation  and  tracking  error. 

Theorem  1  The  velocity  observer  of  (37),  (39),  and  the  control  law  of  (70),  (81)  ensure  that  the 
tracking  error  is  semi-globally  uniformly  ultimately  bounded  (SGUUB)  as  shown  below 


% 


11^)11  <  \/-HWo)ll  exp(-¥fO  + 


\Ae20 


A; 


A-4  A2A3A5 


(1—  exp(— m)) 


(88) 


where 


x 


TlT 


(89) 


and  e0,  A2,  A5  are  positive  constants  and  A3,  A4  are  positive  constants  given  by  the  following  form 


A3  =  nrin{l,mi,  A)02}, 
A4  =  max{l,m2,/c02}, 


(90) 


under  the  condition  that 


A9 

kr  >  3,  kz>  3,  a  >  —to,  k02  >  2 £c4f3, 

Zi 


ko3  > 


m  1 


£c0  +  2£c4  +  Ccl(£2  +  £7)  +  5£cl£8W-p  ||'7(0)||2  +  4  ° 


a3 


A2A3A5 


(91) 

(92) 


where  £c0,  £cl,  £c4  and  e4  to  £3  are  some  positive  constants.  The  details  of  subsequent  stability 
analysis  is  proved  in  Appendix  D. 


Remark  1  According  to  Theorem  1  and  its  subsequent  composite  stability  analysis  from  (116)  to 
(129),  V(rj(t))  is  bounded.  From  (88),  if  the  observer  and  controller  gains  in  (92)  are  selected  to 
satisfy  (121),  it  is  straightforward  to  see  that  r](t)  G  C^.  Hence,  it  ensures  that  the  ep(t),  r(t), 
z(t ),  s(t),  and  x(t)  are  bounded.  We  know  that  all  desired  position  and  yaw  angle  trajectories  are 
bounded,  and  that  R(q)  and  T(0)  are  bounded,  thus  D(R,T )  and  G*(R)  G  C^.  We  can  now  make 
a  conclusion  that  p(t),ip(t),v(t)  G  C^  owing  to  ep(t),ev(t),eg(t),rp(t)  G  £00,  based  on  the  definition 
of  (15),  (19),  (23),  (25).  Due  to  the  boundedness  ofv(t),  we  observe  that  p(t)  in  the  first  equation 

of  (8)  is  bounded,  x  (t)  is  bounded  in  (fl)  because  s(t),  x(t)  G  and  this  ensures  thatp  (£),  0  (t) 
G  Coo  and  we  know  p  (t)  is  bounded  from  (42).  We  also  know  that  v(t),cu(t)  are  bounded  in  (43), 
and  hence,  ev(t),  rp(t ),  z(t)  G  C 00  from  (69),  (68),  (65),  respectively.  Also,  the  expression  for  the 
observed  velocity  v(t)  is  bounded  from  the  definition  of  v(t)  in  (164).  Therefore,  ev(t),rp(t )  G  C 00, 
so  we  know  the  filtered  position  error  r(t)  is  bounded  in  (58).  Owing  to  v(t),  v(t)  are  bounded,  the 
nonlinearity  of  the  aerodynamic  damping  term,  Ni{v),  N\{v),  Y\(y),  Yi(h),  Yi(-)  are  all  bounded. 
Hence,  we  can  state  that  u\(t)  is  bounded  in  (70),  thus  the  translational  control  input  u\(t )  G  C^ 
in  (32).  Owing  to  the  fact  that  u\(t)  is  bounded,  u(i)  is  bounded  from  (31).  From  (65),  we  know 
that  Oft )  is  bounded,  and  hence,  z(t)  G  C ^  from  (64).  Therefore,  since  v(t),  u(t)  have  been  shown 
to  be  bounded,  we  know  that  the  velocity  output  x  (t)  is  bounded  from  (36),  thus  the  aerodynamic 
damping,  N2{oj),  N2{Cj),  Y2(uj),  12( lj),  andY2(-)  are  bounded.  Hence,  x(t),  x0(t )  are  bounded  from 
(42)  and  (40)  because  of  x  ( t )  G  C OQ.  From  the  equations  of  (57),  (56),  (61)  and  (63),  we  can  see 
that  ep  (f),  ev  ( t ),  eg  ( t ),  and  f  ( t )  are  all  bounded  which  can  be  used  to  show  that  uA  (t)  G  £00  * n 
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(79),  thus  Y393  is  bounded  in  (78),  and  consequently  the  torque  input  u2(f)  is  bounded  from  (81) 
using  (64),  (78),  and  (58).  Thus  v  (t),  Cj  ( t )  G  Coo,  in  (12),  G  Coo  from  (194),  and  y(t )  is 

bounded  in  (39).  Therefore  we  can  conclude  that  all  signals  remain  bounded  in  the  velocity  observer 
and  the  closed-loop  system. 

4.1  Observer  Stability  Analysis 

In  order  to  analyze  the  results  of  the  observer,  a  non- negative  function  \f(t)  G  P1  is  defined  as 
follows 

Vo  =  ^ stM*(R ,  T)s  +  ^ k02xTx .  (93) 

The  time  derivative  of  V0  (f)  is 

Vo  =  \sTjf  (. M*(R ,  T))  s  +  stM*(R,  T)s  +  k02xT  £  .  (94) 

Substituting  (54)  into  the  second  term  of  (94)  and  then  arranging  the  equation  yields 

V0  =  \sT  (j2(M*(R,T))+2N*{R,T,x)'j  s  +  sTNf(R,T,xa)s  +  stY*(R,T)s 

- k03sTM *  (R,  T)s  —  k02pxTx.  (95) 

PI  is  used  to  show  that  the  first  term  in  (95)  is  zero  because  of  its  skew-symmetric  property  which 
is  validated  in  Appendix  B,  and  hence,  we  have 

VQ  =  —k03sTM*(R,  T)s  -  k02pxTx  +  stN*(R,  T,  xa)s  +  sTY*(R,  T)s.  (96) 

The  term  N*(R,T,x0 )  is  a  Centrifugal/Coriolis  force  matrix  which  is  derived  from  the  difference 
between  N*(R,T,x)  in  the  modeling  equation  of  (12)  and  N*(R,T,x0)  in  the  auxiliary  signal  y(t) 
in  (39)  of  the  observer  design  as  shown  in  (218)  in  the  Appendix  F.  Based  on  the  structure  of 
N*(R ,  T.  x0 )  and  Al  the  term  N*(R,  T,  xc )  can  be  shown  that 

\\N((R,T,x0)\\  <  £cl  ||i0||  where  ^cl  G  3?1  is  a  positive  constant.  (97) 

According  to  P2,  (14)  in  A3,  and  (97),  (96)  can  be  upper  bounded  as  follows 

Vo  <  -k03mi  ||s||2  -  k02P  \\x\\2  +  £ci  ||i0||  •  ||s||2  +  ^c0  ||s||2  .  (98) 

From  (40),  x0(t)  can  be  upper  bounded  as 

||^o||  <  x  -\-  (3  ||x|| ;  (99) 

therefore,  (99)  can  be  used  to  upper  bound  (98)  uniformly  as  follows 

Vo  <  —(ko3mi  —  ^c0)  ||s||2  —  k02(3  ||5i||2  +  ^ci(  x  +  (3  p||)  ||s||2  .  (100) 
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4.2  Controller  Stability  Analysis 

In  order  to  prove  the  tracking  result,  a  non- negative  function  V\  (t)  6  ft1  is  defined  as  follows 


Vi  (t)  =  +  \rTr  +  \zTjz- 

After  taking  the  time  derivative  of  (101),  substituting  from  (20),  (71),  and  (87),  we  have 

Vi  =  ej  ep  +  rTr  +  zTJz 

=  el  -S(uj)ep  +  —  (rp  -  aep  -  6)  +  (—  -  l)RTpd 
F  m  m 


(101) 


+rT  —krr  +  kr 


- S(u)rp- % 
0 


+zT  —kzz  +  kzuo  +  Y393  —  B^r  +  mS(8)v 


av  +  Yi#i 
0 


(102) 


where  ev(t)  =  rp(t )  —  aep  —  6  was  used  for  ep(t).  After  collecting  terms,  deleting  the  zero  terms, 
canceling  the  cross-terms,  substituting  rp(t)  from  (68),  and  using  the  relationship  given  by 


-S(S)z 

Tz(e)z 


-S(S) 

TAB) 


z  =  B z, 


(103) 


the  equation  ( 102)  becomes 


a  1  1 

Vi  =  —krrTr  —  kzzT  z  —  —e„ep - e„8  +  eE( - 1  )RTpd  +  krr^(mv) 

m  F  m  F  F  m  F 

+rp(av  +  YiOi)  +  kzzTCu  +  zTmS(S)v  +  zTY393. 

After  upper  bounding  the  first  three  terms  in  (104)  and  arranging,  we  have 


V\  <  —kr  Mr* II  —  kz  \\z\ 


+rj  (krm  +  a)v +  Yi9i  +  zT  Y393  +  kzu  +  mS(8)v 


ep\\2  +  ep  (--i )RTpd--8 


(104) 


(105) 


The  first  bracket  term  in  (105)  can  be  upper  bounded  using  the  inequality 

ep  (— -1  )RTpd-—8  <  ||ep||e0  <  \  (^2  ||ep||2  +  (106) 

^  to  to  J  2  \  A  2  / 

An  upper  bound  for  the  second  and  third  bracketed  terms  in  (105)  is  now  sought.  The  definition 
of  $(fi,d))  in  (43)  is  upper  bounded  as 

$(v,u>)  =  D(R,T)  x<  \\D(R,T)\\  x  =  d0  x  (107) 

where  do  G  'R1  is  a  positive  constant.  Hence 

||fi||  <  do  x  and  ||cD||  <  do  x  .  (108) 

Yi(v)  and  Y3(v,uj)  are  bounded  utilizing  (72)  and  A3  as  follows 

Yi9i  <  £c2  x  and  Y393  <  £c3  x  (109) 
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where  £c2,  £c3,  are  positive  constants.  From  the  definition  of  s(t)  in  (41),  x  (t)  can  be  upper  bounded 
as 


x 


x 


<  INI+/?||NI 

<  2  ||s||2  +  2/32  ||x||2 . 


(110) 

(111) 


Upper  bounds  for  the  last  two  bracket  terms  in  (105)  can  be  expressed  in  the  following  form  by 
utilizing  (108)  for  v(t)  and  L d(t),  (109)  for  Yi91  and  Y393  ,  (110)  for  x  ,  and  (111)  for  ||x||  to 
produce 


Vi(t)  <  -kr\\r\\2  -  kz\\z\\2  -  ( 


1 


2  j)  IM2  + 


m 


+3||r||2+^2 

X 

2 

+  a 2 

X 

2 

+3|kl|2  +  e23 

X 

+  k2z 

X 

+  ( krm ) 


x 


m2\\sm' 


x 


After  rearranging  (112),  we  have 


(112) 


Vi  <-{kT-  3)  j | r 1 1 2  -  (kz  -  3)  ||2||2  -  (—  -  h)  j|e  f  +  (  4  ||i 

m  2 


±  <113> 


where 


I2 

Uoo 


icA  —  Cc2  +  £c3  +  «2  +  ( krm )2  +  k~  +  m2  ||5'(5) 

After  using  (111),  (113)  becomes 

U  <-(kr-  3)  ||rf  -  (*,  -  3)  ||2||2  -  (-  -  h)  ||Cp|P  +  2ecl(||s||2  +  f  ||i||2)  +  A. 

Ill  Z  Z  A  2 


(114) 


4.3  Composite  Stability  Analysis 

The  performance  of  the  proposed  controller  and  observer  can  now  be  examined  by  combining  the 
non-negative  functions  V0(t)  in  (93)  and  V± (t)  in  (101)  as  follows 

V  =  V0  +  Vi  =  l-eTpep  +  Vr  +  l-zTJz  +  l-sTM*(R,  T)s  +  h02xTx.  (115) 

The  composite  function  V(t)  now  has  the  property 

Ja3  ||r?||2  <  V  <  Ja4  \\r]\\2  (116) 


where  (89)  is  used.  After  taking  time  derivative  of  (115),  and  utilizing  the  bounds  on  V0  (t)  and 
Vi  (t)  from  (100)  and  (114),  we  have 


V 


~(kr  -  3)  ||r||2  -  (kg  -  3)  |N|2  - 


—kwfl  ||NI  +  £, 


cl 


X 


+  P\\x 


(£  -  y)  llepl|2  -  (kosTTH  -  £c0)  INI2 

ill  z 

llsl|2  +  2^c4  INI2  +  2^c4/32  |N||2  + 


(117) 
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We  can  now  upper  bound  V{t)  of  (117)  as  follows 


v  <  -(fc-3)||r||2-(fc«-3)||2||2-(^-^)|M: 

III  Z 


^03™1  £co  Cel 
+  (2^c4^2  —  ^02/5)  p||2 


X 


^11*11-2^  m 


2Ao 


(118) 


X 


in  (118)  can  be  upper  bounded  by  using  (190)  from  Appendix  D  and  substituting  (110)  into 


(190)  as  follows 


<  (2  +  £3)  HI  +  (2  +  £3)/?  ||x||  +  £4  ||z||  +  (£1  +  £5)  ||r||  +  (ck£i  +  £g)  ||ep||  +  £2  +  £7- 


(119) 


After  substituting  (119)  into  (118)  and  arranging,  we  have 


V  <  -{kr  -  3)  ||r||2  -  (kz  -  3)  ||zf  -  (—  -  ^)  ||ep||2  -  [fc03mi  -  £c0  -  2£c4]  IM|2 

Tfl  Z 

—• Cd  ((2  +  £3)  INI  +  (3  +  £3 )p  ||x||  +  £4  ||z||  +  (£1  +  £5)  ||r||  +  (a£i  +  £6)  ||ep||  +  £2  +  £7)  ||'sl|" 
+mc4p2-k02p)\\x\\2  +  ^-.  (120) 

Utilizing  (89)  for  the  second  line  in  (120),  V(t)  yields 

V  <  -(kr-3)\\r\\2  -  (kz-3)\\z\\2  -  -^)\\ep\\2 

III  Z 

-  [kosrri!  -  £c0  -  2£c4  -  £cl(£2  +  £7)  -  5£cl£8  |M|]  ||s||2 

~(k02p  -  \\x\f  +  (121) 

where 

£8  =  max{(2  +  £3),  (3  +  £3)/?,  £4,  (£1  +  £5),  (ck£i  +  £6)},  (122) 

An  upper  bound  can  be  written  for  (121)  as 

V  ^  -A5  [||r||2+ ||^||2  +  ||ep||2  +  ||s||2  +  ||5||2]  +  ^ 

<  -A5|H|2  +  ^  (123) 

where  a  positive  constant  scalar  A5  G  ’ft1  is  given  by 

A5  =  min  {{kr  —  3),  (kz  —  3),  (^  —  4f),  (k02P  —  2£c4/d2), 

[k03mi  —  £c0  —  2£c4  —  £cl(£2  +  £7  +  5£8  llhll)]}  •  (124) 


provided 

ko3  >  ^  [£co  +  2£c4  +  ^ci(£2  +  £7)  +  5£c1£8  \\r]\W  ,  (125) 

and  the  conditions  for  gains  in  (90)  and  (91)  are  met.  A  sufficient  condition  for  (123)  is  found  from 
(116)  using 


to  write 


(127) 


V<-2^V+it, 


2A2  : 


provided 


ko3  > 


m  i 


£c0  +  2£c4  +  Ccl(^2  +  £7)  +  5^c1^8\/ (t) 


(128) 


which  yields  a  new  sufficient  condition  for  (123).  Solving  the  differential  inequality  in  (127)  yields 


r  <  l/(0)exp(-^f)  +  (1  -  exp(-fy))  . 


Then  we  can  write 

and  from  (116)  we  can  write 

and  combining  these  yields 


V  <  1/(0) 


r(0)  <  f  ||,(0) 


A4^q 
4A2A5  5 

2 


(129) 

(130) 


V  <% 


Mo)  I 


A4£q 

4A2A5 ’ 


(131) 


which  can  be  combined  with  (128)  to  produce  the  sufficient  condition  for  given  in  Theorem  1  to 
satisfy  (123).  Substituting  (129)  into  (126)  yields 


INI  <  Vfe  IM0)l|Mxp(-^)  +  ^(l-exP(-31p)). 

Therefore,  the  result  of  Theorem  1  can  be  obtained. 


(132) 


5  Simulation 


The  output  feedback  tracking  control  was  simulated  using  a  small  quad-rotor  unmanned  aerial 
vehicle  [13]  as  depicted  in  Figure  1.  The  inertial  parameter  of  the  simulation  vehicle  are  borrowed 
from  [14] 


m  =  0.9  [kg\,  J  = 


0.32  0  0 

0  0.42  0 

0  0  0.63 


[kg  ■  m2 


where  g  =  9.81  [to/  sec2]  is  acceleration  of  gravity,  and  are  assumed  to  be  constant  while  following 
the  trajectory.  The  desired  position  and  yaw  trajectory  are  given  in  the  following  form,  respectively 


Pdx 

Axsin(wt)(l  —  e(  B^) 

Pd(t)  = 

Pdy 

= 

Ay  cos(wt)(l  —  e<'~Byt3'1) 

.  Pdz  . 

Az(  1  - 

[tpd]  =  [  c0sin(2^t)  ]  (rad) 

where  Ax  =  Ay  =  1,  Az  =  1,  Bx  =  By  =  Bz  =  5,  T  =  2.5,  w  =  2  •  2^, and  cQ 
position  and  angle  of  the  quad-rotor  at  the  center  of  mass  are  selected  as  follows 


(133) 

(134) 
1 .  The  initial 


p(0)  =  [0.1,0.1,0.1]T,  0(0)  =  [0.1,  —0.1, 0.1]r,  x(0) 


[0,0,  -0.1,0,0,0]T 


where  all  parameter  estimates  x(t)  are  initialized  to  zero  except  the  position  of  z-direction  (j3z(0)  = 
—0.1).  The  initial  orientation  and  constant  vector  are  chosen  as  follows 

^(O)  =  [1, 0, 0, 0]T,  <5  =  [0,0,<53]t,  S3  =  -l. 
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The  constant  control  parameters  for  observer  and  controller  were  iteratively  chosen  to  be 

&01  =  4,  ko 2  =  4,  /c03  =  2,  (3  =  2,  /c03  =  4, 

=  5,  /cr2  =  1,  kz  =  45,  a  =  50, 

which  satisfies  the  condition  of  (92). 

Figure  3  shows  the  position  tracking  of  the  quad-rotor  to  the  desired  trajectory  Pd(t).  The  actual 
quad-rotor  trajectory  represented  by  the  solid  line  follows  the  desired  trajectory  represented  by  the 
dotted  line  which  is  commanded  to  go  up  1  [m]  high  and  rotate  around  circular  orbit  of  radius  of  1  [m] 
in  the  plane.  Figure  4  shows  position  tracking  at  each  axis  corresponding  to  the  motion  in  Figure 
3  and  the  last  one  shows  yaw  tracking  result  to  the  desired  trajectory  Figure  5  represents 

the  position  errors  about  the  each  coordinates  ( x,y,z )  and  yaw  angle  errors.  Figure  6  shows  the 
control  inputs.  The  translational  force  input  u\  (t)  is  collectively  steady  when  the  UAV  rotates  at 
the  orbital  trajectory.  The  torque  commands  u^ft)  periodically  changed  when  they  rotate  around 
the  circle.  Figure  8  shows  the  estimated  output  of  velocity  observer. 


6  Conclusion 

The  goal  of  designing  and  output  feedback(OFB)  controller  for  a  quad-rotor  UAV  system  has  been 
demonstrated  mathematically  and  via  a  computer  simulation.  The  mathematical  result  shows  that 
a  semi-global  uniformly  ultimate  bounded  (SGUUB)  tracking  result  is  achieved.  The  nonlinearities 
of  the  damping  term  were  included  in  the  system  modeling  and  it  was  linearly  parameterized 
because  the  velocities  or  other  factors  are  assumed  to  be  unmeasurable.  While  the  output  feedback 
control  design  was  predicated  on  a  hypothetical  sensing  system  that  only  produces  angular  and 
linear  positions,  it  does  appear  that  low-cost  GPS  or  a  camera  based  units  may  provide  justification 
of  this  approach.  It  worth  noting  that  the  output  feedback  design  has  its  advantage  over  the  full 
state  feedback  where  a  full-state  feedback  controller  could  be  considered  a  less  complicated  subset 
of  the  current  work.  We  believe  this  to  be  the  first  paper  to  present  such  a  comprehensive  result  for 
quad-rotor  tracking  control  based  on  only  position  measurements.  Although  the  focus  of  this  work 
is  the  quad-rotor  class  of  aircraft,  the  results  are  directly  applicable  to  other  aerial  vehicles  such  as 
the  co-axial  helicopter  and  to  unmanned  underwater  vehicles. 
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A  Development  of  Dynamic  Model  in  the  Inertial  frame 

The  dynamic  equation  of  (1)  describing  the  dynamics  of  the  quad-rotor  can  be  directly  written  in 
the  form 

Md>(7,  V)  =  C(R,  T)D(R,  T)x  +  h(R,  T,  x)  +  G{R)  +  BU  (135) 


23 


using  the  definition  of  $(n,o;)  in  (10)  and  the  following  substitutions 


M  =  \  °3f  ]  e  5£6x6,  C(R,T)  =  \  °3f,  ]  e  d£6x6, 

1^3x3  J  J  [  ^3x3  b(JuJ) 

G(R )  =  [  Ci\K]  }  6  S6,  h(R,T,i)  =  [  1  6  S6,  (136) 

<-/3xl  J  |_  iV2la;J 

B  =  \  °3x3  ]  e  d£6x4,  U  =  \  Ul  ]  e  3?4  . 

O 3x1  B 2  J  L  . 

The  time  derivative  of  $(n,o;)  in  (10)  can  be  related  to  x(t)  by  differentiating  (9)  and  applying  the 
definition 

^  (D(R,T))  =  D(R,T,x)  (137) 

where 

jtWR’T))  =  [*<£?  s(i?T)  =  J?T  =  -SMflT“d 

jt{T~ 4(0))  =  A(r-1(0))0  e  3^3x3  where  ^(T-l(0))  e  ^3x3x3  is  a  tensol.  (138) 

to  yield 

4>(h,  u)  =  y  (D(R,  T))x  +  D(R,  T)x  =  D(R,  T,  x)x  +  D(R,  T)x.  (139) 

LLL 

Multiplying  by  M(-),  substituting  from  (135)  for  MQ(y,cu)  and  from  (9)  for  $(n,a;),  and  arranging 
terms  yields 

MD(R,  T)x  =  C(R,  T)D(R,  T)x  -  MD(R,  T,  x)x  +  h(R,  T,  x)  +  G(R)  +  BU.  (140) 

It  can  be  shown  that  Dt(R,T )  =  [R(q),  03x3;  03x3,  T~t^  <E  di6x6.  After  multiplying  (140)  by 
Dt  =  Dt(R,  T ),  we  have 

DtMD(R,  T)x  =  Dt  ( C(R ,  T)D(R,  T )  -  MD{R ,  T,  x))  x  +  DTh(R,  T,  x)  +  DrG(R )  +  DTBU. 

(141) 

Equation  (141)  is  now  written  in  the  compact  form 

M*(R,  T)x  =  N*(R,  T,  x)x  +  h*(R,  T,  x)  +  G*(R)  +  B*(R,  T)U  (142) 

where  the  corresponding  matrices  were  substituted  as  follows 


M*(R,  T)  4  Dt(R,T)MD(R,T)  =  [mh,03x3;03x3,T-T(e)JT-\e)}, 

N*(R,T,x)  =  Dt(R,T)  (C(R,T)D(R,T)  -MD(R,T,x))  in  (209) 
h*(R,T,x)  =  DT(R,T)h(R,T,x)  =  [RN1(v)]T-t(Q)N2(uj)},  (143) 

G*(R)  =  DT  (R,T)G(R)  =  [RG(R);Osxi], 

B*(R,T)  =  Dt(R,T)B  =  [RB1,03x3-,03x1,T-t(Q)B2]. 


24 


B  Proof  of  Model  Property  PI 


The  validity  of  the  skewed  symmetric  relationship  in  PI  is  shown  below.  The  definitions  of  M*(R,  T ) 
and  N*(R.  T.  x)  from  (143)  in  Appendix  A  are  first  applied  to  the  matrix  term  of  PI,  that  is 

(. M*(R ,  T))  +  2N*(R,  T,  x) =  0,  V  £  e  5i6 


to  yield 

(|  (. M*(R ,  T))  +  2 N*(R,  T,  x))  £  =  f  (|  T)MD(R,  T)))  f  (144) 

+eT  (2ZlT(i?,  T)  (C(i2,  T)D(R,  T)  -  MD(R,  T,  x))) 

The  definition  of  Z)(i?,  T,  x)  from  (137)  in  Appendix  A  is  now  applied  and  terms  collected  to  yield 


£t  (M*(iZ,  T))  +  2N*(R,  T,  x))  £  =  £T  (2HT(i?,  r)C(i2,  T)D(R,  T))  £  (145) 

+£T  ( DT(R ,  T,  x)MD(R,  T)  -  DTMD(R,  T,  x))  £ 

where  M(-)  is  a  constant  matrix  and  hence  =  0  has  been  used.  It  is  now  useful  to  invoke  the 
symmetric  property  of  M(-)  to  write 


e  (. DT(R ,  T,  x)MD(R,  T))  e  -  e  (- DTMD{R,  T,  x))  £  =  0, 


which  allows  (145)  to  be  written  as 


(146) 


e  (M*(R,  T))  +  2N*(R,  T,  x)  j  £  =  2?  DT(R ,  T)C(R,  T)D(R,  T)£.  (147) 

It  is  now  possible  to  introduce  a  new  vector  e  TG  defined  as  =  D(R.  T)£  in  order  to  rewrite 
(147)  as 

(  4  T))  +  2N*(R>  T>  C  =  2 C'TC(R,  T)?  (148) 


f  ^-(M*(R,T))  +  2N*(R,T,x)jt  =  2ZuC(R,T)S'  (148) 

where  it  is  clear  that  if  C(R,T)  is  skew  symmetric  then  the  right-hand  side  of  (148)  is  zero.  The 
definition  of  C(R,T)  in  (136)  is  substituted,  and  partitioned  into  subvectors  to  write  (148)  as 


fC(R,T)(=[g  ftT]  ofiT*  S<X>)  I  =  s(“)&  +  &  (149) 

where  it  is  clear  that  the  skew  symmetry  property  of  S(u)  and  S(Juj)  can  be  invoked  to  write 


and  hence, 


2i,TC(R,T)il  =  0, 


y  -  (M*  (B,  T))  +  2iV*  (B,  T,  i)  U  =  0 


(150) 


(151) 
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C  Development  and  Derivative  of  Control  Signal  u\(t) 

The  term  of  interest  is  first  repeated  from  (30)  and  the  definition  of  f  ( t )  in  (27)  is  inserted  as  a 
reminder  to  form 


—S(u)rp  +  ^ ev  —  RT  Pd  +ct(^  —  1  )RTPd  +  G(R) 

4- 

'  N^v)  ' 

4- 

—S(8)u  +  BiUi 

—uzd 

0 

r2(©v 

~  ^  .  (152) 

With  Bi  =  [0,0, 1]T,  it  can  be  seen  that  the  control  input  U\  (t)  will  only  act  on  the  z-axis  trans¬ 
lational  dynamics.  A  new  control  signal  will  be  injected  that  creates  control  signals  to  all  three 
translational  axes;  these  injected  control  signals  then  become  the  tracking  objectives  for  the  rota¬ 
tional  dynamics.  The  control  term  from  (152)  can  be  divided  into  two  parts;  B M(-)  defined  in  (34) 
and  flit)  G  5i4,  and  hence,  becomes  as  follows 


— S(6)uj  +  BiU\ 

i 

Co 

1 

1 _ 

c 0 

Tz(Q)u 

_  TZ(Q)  0  J 

_  Ul  _ 

(153) 


—T  - 

The  term  ui(t)  =  [ul  ,ui]  G  d£4  with  u\  ( t )  G  !R3  can  be  added  and  subtracted  to  /i(t)  of  (153)  as 
follows 

r  .  i  r  .  i  r  -  i  -  i 

(154) 


UJ 

UJ 

Ul 

UJ—  Ul 

-  Ul  . 

=  Ui  + 

-  U1  _ 

-  W  . 

=  Ul  + 

0 

p  =  —  f/  i  -  Ui  + 

The  term  z(t )  can  be  defined  to  be  a  measure  of  the  closeness  of  uj(t)  to  u\  ( t )  as 


z  =  UJ —  u  i, 


which  can  then  be  written  as 


z  =  oo  —  Bzu  i  where  Bz 


10  0  0 
0  10  0 
0  0  10 


(155) 


It  is  interesting  to  note  that  if  z{t )  is  zero  then  the  control  signal  Ui(t)  has  been  effectively  injected 
into  the  open- loop  filtered  tracking  error  dynamics.  The  manipulation  of  //(f)  can  be  continued 
from  (154)  using  the  definition  of  z(t )  to  yield 


A i  =  «1  4-  Q  • 

After  multiplying  / x(t )  of  (156)  by  B^{-)  of  (153),  we  have 


B^l  =  B^ui  +  B M 


(156) 


(157) 


and  thus 


—S(6)u  +  B\U\ 
Tz(6)w 


B^U\  +  B n 


z 

0 


(158) 


The  derivative  of  a  \  (t)  is  obtained  by  taking  time  derivative  of  (66).  It  is  perhaps  more  illustrative 
to  abbreviate  (66)  as 

hi  =  B^U  (159) 
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where  B  1  ( 


is 


B"  = 


and  hence,  after  differentiation,  we  have 


0  -f  0  0 

k  000 

_ LM  n  n  — 

S3  c  <t>  c  </> 

0  0  10 


(160) 


(161) 


where  is  now  clear  that  two  derivatives  are  required.  Term  by  term  differentiation  of  B  1  (•)  yields 


d 


(b:1)  = 


dty  * 


0  0  0 
0  0  0 


0 

0 


__  1  0  n  n  -c(frs0’d+s(frc0(j) 

6 3  c 2(j)  c2</> 


0  00 


0 


(162) 


where  0(f)  is  measurable  and  <j)  (t)  =  Tx(Q)u  and  9  (t)  =  Ty(Q)Cj  are  estimated  from  (35).  Next, 
the  time  derivative  of  U(t)  is  rewritten  from  (80)  as  follows 


drr  ,  r 

—U  =  —kr  r  + 
dt 


d_ 

dt 


1 


—  ev  —  —  e„ 

m  u  m  B 


0 


S{u>)  (%RTPd  ~  uRTpd  -  RT  Pd  +mgRTEz  j  +  a  (l  -  £)  RT  Pd  +RT  Pd 

dJzd 

where  (56),  (57),  and  (63)  are  utilized.  The  term  Yi(y)9\  is  a  general  representation  of  the  nonlinear 
aerodynamic  damping,  and  therefore,  the  time  derivative  cannot  be  written  explicitly  until  a  specific 
model  has  been  assumed  (such  a  model  is  assumed  for  use  in  the  Simulation  section  and  the  resulting 
derivative  ^  (Yi(f))0i)  is  calculated  in  Appendix  E. 3.  Therefore,  U\  (t)  is  obtained  by  substituting 
(162)  and  (80)  into  (79). 


D  Details  of  Stability  Analysis 


Development  of  an  upper  bound 


x 


in  (118)  requires  bounds  for  p  (t)  and  0  (t). 


D.l  Upper  Bound  for  P  (t) 

As  a  starting  point  for  the  bound  on  p  (t),  ev(t)  in  (19)  is  substituted  into  the  definition  of  rp(t)  in 
(25)  and  the  result  solved  for  v(t)  to  yield 

v  =  —  (rp  -  aep  -  6  +  RTpd)-  (163) 

m 

We  can  now  use  (163)  and  v(t)  =  v  —  v  to  obtain 

v  =  -v  +  —  (rp  -  aep  —  6  +  RTpd )  .  (164) 
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Now  p  (t)  can  be  expressed  by  using  (35)  and  substituting  from  (164),  as  follows 

p  =  Rv 

=  —Rv  H - R  (rp  —  aep  —  b  +  RTpd )  , 

and  then  we  utilize  (35),  (165),  and  p  ( t )  =  Rv  to  yield 

p=  —  p  +—R  (r„  -  aep  -  6  +  RTpd )  ■ 
m 


(165) 


(166) 


The  triangle  inequality  can  be  utilized  to  create  an  upper  bound  for  p  (t)  in  (166)  as  follows 


< 


P 


R 

m 


\rp\\  +  a 


R 

m 


Ml  + 


R 

m 


||<5  +  RTpd\\ , 


which  can  be  further  bounded  as 


< 


P 


+  £i  pp  +  a£\  e„  +  £2 


(167) 


(168) 


where 


£1  -  -J-  SUP  II^IUoo  ’  £2  -  “  SUP  H^IUoo  ( ll^ll  +  sup  \\RT\L  sup  IMII )  •  (169) 

m  ye  m  ye  \  ye  \/t 


D.2  Upper  Bound  for  ©  ( t ) 

To  begin  the  development  of  a  bound  for  0  (t),  £j(t)  =  uj  —  Cj  is  solved  for  c v(t)  producing 

Cj  =  uj  —  6j. 

After  multiplying  (170)  by  the  transformation  matrix  T ( 0 ) ,  we  have 

T(0)  u  =  T(Q)u-T(Q)u, 

which  is  equivalent  to 

0=  T(0)u>—  0  . 

The  definition  of  z(t)  in  (31)  is  solved  for  u>(t)  to  yield 

c u  =  z  +  Bzui , 

and  the  control  ui(t)  in  (66)  is  substituted  into  (173)  to  yield 


(170) 

(171) 

(172) 

(173) 


lv  =  z  +  BZB.,  —krr  +  k. 


mv 


+ 


^  nr 1  1  n  r  |  g 

■(rp  —  aep  —  6)  +  av 


m  ' 


—  —  Gt 
m  \  m 


(^  -  1  )RTpd  +  RT  Pd  —G(R)  -  17(7)01 

Uzd 


(174) 
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where  rp(t)  =  mv  was  used  for  r(t),  ev(t)  was  substituted  from  (25),  and  rp(t)  was  substituted  from 
(69).  By  substituting  (174)  into  (172),  0  (t)  can  be  expressed  as  follows 


0  =  -  0  +Tz  +  TBzB~l  I  - krr  +  k 


/j, 

■{rp  —  aep  —  6)  +  av 


mv 

0 


„.t  i 
m  "  v  m 


-  1  )RTpd  +  RT  Pd  - G(R )  -  Yi(fl)0i 

Uzd 


(175) 


In  order  to  combine  the  two  v(t)  terms,  the  matrix  Bj(-), 

e  Ji4x3, 

and  the  equality  RT  p  (t)  =  v(t)  from  (35)  are  used  to  formulate  the  equality 


BTZ  = 


h 

Olx3 


(176) 


v 

0 


=  Blv  =  BjRTp, 


(177) 


which  is  then  substituted  into  (175)  to  yield 


0=  -  0  +Tz  +  TBzB~^Btz  ( krm  +  a)RT  p  -TBzB~lkrr 

-r  I  ~%rP  -  Y1(v)61  +  ^(«2  -  1  )  +  %6-  a(£  -  l)RTpd  +  RT  Pd  -G(R) 

Uzd, 


+TBZB, 


.(178) 


From  the  A2,  the  nonlinear  aerodynamic  damping  term  Yj  {v)0\  in  (178)  can  be  bounded  by  the 
following  expression 

lin(h)0i||<  ^li¬ 


lt  can  be  used  to  show  v(t)  =  RT  p  which  can  be  substituted  into  (164)  to  show  that 

~  ~  1  1  /  c  j.'!  .  \ 

v  =  —R  p  H — rp - ep - (o  -  R  pd). 

m  m  m 

v(t)  can  now  be  upper  bounded  in  the  same  manner  as  ||p||  in  (167)  as  follows 

1  „  „  a  „  „  1 


(179) 


(180) 


7  <  R1 


P 


~  Vp  +-  eP  +- 
to  to  to 


S  +  RTpd  . 


The  definition  of  r(t)  in  (26)  leads  to  the  bound  on  rp(t)  given  by 


I  rp  <  |  r 


The  bound  in  (179)  can  now  be  upper  bounded  in  the  manner 


in  Will  <e 


c4 


RJ 


p 


—  IWI  +  —  IMI  +  —  ||<5  +  RTpd\ 

TO  TO  TO 


(181) 


(182) 


(183) 
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by  using  (181).  An  upper  bound  can  now  be  developed  for  0  (t)  from  (178)  as  follows 


0 

< 

0 

T ||  •  \\z\\  +  [  TBzB~'Blz  ( krm  +  ajR1  +  TBzB~'Blz  •  £c4  R 1 


tbzb„'bT(—  +  k 


a 


tbzb~xbt^—  -  — ; 

M  K  m  m 


m 

o ?  1 


I  TB,BZxBtM  ^ 
m 


(y 

-U  \\tb-.b-'b 

to  p 


-1  uT 

z 


\\TB,B-'Bl 


cv  1 

-6  +  o:(  -  l)i?Tpd  +  i?T  Pd  +G(R) 

m  m 


V 


(184) 


+  \\TBzB^B3\\  ■  \\uzd\\  +  \\TBzB~1Btz\\  •  ^  \\6  +  RTpd || 

where  the  abbreviation  T  =  T(0)  and  B3  =  [0, 0,  0,  l]7  for  matrix  dimension  were  used,  and  (182), 
(183)  have  been  utilized.  The  following  constants  are  introduced  to  simplify  expression  (184) 


£3 

£4 


=  £4£g  BZB  Bz  (krm  +  a)  +  £4£9  BZB  Bz  .  £ 


5  c4 


sup  iiml, 

ye 


£5  —  £4 

£q  =  £4 


a 


BzB~lBz  (  b  k, 


m 

,a 2  1 


+  ^£4i|5ZJB-iJEqi. 

/m  II  z  1 1  too 

zoo  in 


BZB~  Bi - ) 

M  v  m  m 


a 


-Uet  \\BzB^Bt' 


m 


£7  — 


£4  BzB~lBTz  . 

A1  z  *oo 


a 


m 


I  <5 II  +  ck 


1 


-  1 


TO 


M  ^  II ioo 


£9  Slip  1 1  1 1  +  £9  Slip 

Vi  Vi 


Pd 


+£4  ||TOr?,71S3||ioo  sup  |K,||  +  (^£4 


inT 

M  ^  ioo 


(185) 


1*11  +  £9  s«p  llpsl 

it 


£9  — 


sup  ir 

ye 


thereby,  creating  the  final  bound  for  0  (t)  as  follows 


• 

0 

< 

0 

+  £4 IM  +  £3 

p 

+  £5  Ml  +  £6  kp||  +  £7- 


(186) 


The  definition  of  x  (t)  is  formed  from  the  right-hand  terms  in  (35)  and  (36)  and  an  upper  bound 
applied  as  follows 


X 

_ 

V 

< 

p 

+ 

0 

.  0 

(187) 


The  bound  in  (168)  for  p  (t)  and  the  bound  in  (186)  for  0  (t)  can  now  be  substituted  into  (187)  to 
create  an  upper  bound  for  x  (t)  in  the  manner 


(189) 


Additionally,  the  terms  p  (t)  and  0  (t)  can  be  individually  bounded  as 


p 

< 

X 

0 

< 

X 

Now  the  substitution  of  the  bounds  from  (189)  and  (182)  into  (188)  yields 


x 


<(2  +  £3) 


+  £4  II^H  +  (^1  +  £5)  IMI  +  (o£  1  +  £q)  \\ep\\  +  £2  +  £7- 


(190) 


E  Simulation  Notes 

E.l  Desired  Trajectory 

The  desired  trajectory  along  the  three  linear  directions  is  given  as  follows.  The  desired  trajectories 
for  the  :r-axis  are 

pdx  =  Ax  sin(wt)(l  —  e*-0'5^) 

pdx  =  Ax(w  cos(wt)(l  —  +  sm(wt)l.bt2e^°'5t^) 

Pdx  =  Ax(-(w)2  sin(wt)(l  -  e(-aM3))  +  2wcos(wt)lM2e^0M^ 

+  sin(wt)3te(-aM3)  -  sin(wf)(1.5f2)2e(-a5t3)) 

Pdx  =  Ax(— (w)3  cos(wf)(l  —  e^0"54'3))  —  (w)2  sin(wt)1.5t2e('_0"5i3')  (191) 

— 2((ry)2sin(wt)1.5t2e('_0'M3')  +  3w  cos(wt)te('~0M^ 

—w  cos(wt)(1.5t2)2e('~0M3^)  +  (3w  cos(wt)te('~0M3^  +  sin(zet)3e('^0'M3) 

— 4.5  sin(wf)f3e^a5t3))  —  rocos(ryt)(1.5t2)2e(-_0'M3->  —  2  sin(zcf)(1.5f2) 

3te0O.5t3)+sin(wt)(L5t2)3e(-O.M3)); 

for  the  y-axis  are 

Pdy  =  Ay  cos(wf)(l  —  e(— °-5*3)) 

Pdy  =  Ay(—wsm(wt)(l  —  e('~05t3^)  +  cos(wt)lM2e('~05t3^) 

Pdy  =  Ay(—(w)2  cos(wt)(l  —  e^0M^)  —  2wsm(wt)1.5t2e('~0M3^ 

+  cos(wt)3te^0M^  —  cos(rcf)(1.5f2)2e^°'5t3)) 

Pdy  =  Ay((w)3  sin(ret)(l  —  e(— °-5*3)))  _  (w)2  cos(u'f)1.5f2e<~0,5t3') 

— 2((w)2  cos(ret)1.5t2e('^0"M3')  +wsm(wt)3  te^0’5^ 

— w  sin(rf;t)(1.5t2)2e^_0'5t3->)  —  w  sm(wt)3te^0M^  +  cos(wt)3e('~0M^ 

—  cos(wt)3t(1.5t2)e('^0M 3-)  +  wsin(n;t)(1.5t2)2e^0'M3')  —  cos(uh)2(1.5f2) 

3te(-0.5t3)+cos(wt)(L5t2)3e(-0.5t3))) 

and  for  the  z-axis  are 

pdz  =  Az(l  -  e(“aM3)) 

pdz  =  Az1.5t2e^ 

Pdz  =  Azl.b(2t  -  1.5f4)e(-°-5t3) 

Pdz  =  A21.5((2  —  6f3)  —  (2t  —  1.5t4)1.5t2)e(_0,M3). 
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E.2  Nonlinear  Aerodynamic  Damping  Term 

The  nonlinearities  of  the  aerodynamic  damping  terms  Ni (v),  N2(w)  were  included  in  the  system 
modeling  (1)  and  specific  knowledge  of  the  system  model  is  given  by 


d\  +  d2  |ui  | 

0 

0 

Vi 

Ai(u)  = 

0 

d3  +  d, !  \v2\ 

0 

V2 

0 

0 

ds  +  d§  u3 

.  V3  . 

9i  +  92  |wi| 

0 

0 

u  1 

N2(w)  = 

0 

93  +  94  <^2 

0 

0J2 

0 

0 

95  +  96  <^3 

UJ3 

(192) 


E.3  Linear  Parameterization 


The  assumed  model  for  the  aerodynamic  damping  terms  in  (192)  were  used  to  create  the  parame¬ 
terization  Yi(v)9i,  Y2(uj)92  by  utilizing  A2  and  invoking  the  estimated  velocities  h(f),  w(t )  given 
by 


Ti(fi) 

0i 

Y2(cj) 

02 


Vi  Vi  ■  |  fii|  0  0  0  0 

0  0  h2  v2  ■  |h2|  0  0  , 

_  0  0  0  0  h3  h3  •  |h3|  _ 

[du  d2,  d3,  d4,  d5,  d6]T  =  [0.065,  0.065,  0.065,  0.065,  0.065,  0.065]T,  and  (193) 

d>l  Cj  1  •  |d>i|  0  0  0  0 

0  0  Cj2  Cu2  ■  \Cj2\  0  0  , 

0  0  0  0  ch3  c53  •  |a>3| 

[gi,  92,  <73, 94,  g5,  9s]T  =  [0.065,  0.065,  0.065,  0.065,  0.065,  0.065]T 


where  the  numerical  values  for  the  damping  parameters  are  borrowed  from  [15].  To  show  the  linear 
parmaterization  of  the  equation  (78)  we  first  need  to  complete  the  derivative  of  u\  (t)  in  (79).  The 
derivative  Yi(h)6fi  in  (80)  to  match  the  derivative  of  the  control  input  ui(t),  u\  ( t ),  which  is  used 
for  torque  input  u2{t )  in  (81)  can  be  found  as  shown  below 


hi  (|hi|  +  hiS5m(hi))  hi  0  0  0  0 

h2  (|h2|  +  v2sgn(v2))  h2  0  0 

h3 


0  0 

0  0 

hi  2  |hi|  hi  0 

0  0  h2 

0  0  0  0 

+  2 d2  ■  |hi|)  hi 
(d3  +  2c?4  •  1 7)2 1 )  h2 
(<4  +  2d§  •  |h3|)  h3 

2 d2  •  |hi| 


0  0 

0  0  0 

2  |h2|  h2  0  0 


di 

0 

0 


h3  2  |h3|  h3 


9 1 


0 

d3  +  2c?4  •  |h2 1 
0 


0 

0 

c?5  +  2g?6  •  |  h3 1 


hi 

h2 

- 

1 

CO 

_ 1 

h3|  +  vzsgn(y3))  h3 


9i 


(194) 


=  H,  v 
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—  Hi  (  —S{£j)v  H - ( BiUi  +  G  —  Yi(v)0i) 


m 


—  77i  S(u))v  H - +  G)^)  —  — 1  i(,O)0i 

\  m  )  m 


where  ifi(-)  e  d£3x3  is  defined  as 


#  i  = 


d\  -f*  2d2  ■  |hi|  0  0 

0  G?3  +  2g?4  -  1 4)2 1  0 

0  0  g?5  +  2g?6  •  1^31 


(195) 


and  |h|  =h  sgn(y)  and  the  first  equation  of  the  modeling  equation  were  utilized.  Substituting 
ui(t)  in  (66)  into  r  (t)  in  (63)  yields 

r=  —Ayr  + 


-S(u)r.  1 

0 


P  meP 


+  Bn 


z 

0 


(196) 


Then  we  substitute  for  U(t)  from  (66)  and  for  4:U(t)  from  (80)  and  also  substitute  general  error 
definitions  from  (56),  (57),  (196),  and  the  derivative  of  Yi(v)9i  from  (194)  to  yield 


d 


u,  =  ft(B^)^-krr  + 
~(B~1)kr  (- krf  + 

“(V 
-(V 

-(V 


a 


(i  -  1  )RTPd  -  G(R )  -  Fi(fi)di  +  RT  Pd  -±ep 


-S(u)rp-±ep 


0 


+  Bl 


Uzd, 

Z 

0 


~  y—S(u)ev  +  G(R )  +  Yi{y)9\  —  RT  i>d  +B\U\ 

0 

Hx  (~S(u>)v  +  ±(BlUl  +  G(R)))  -  H1YYl(v)01 

0 

m  (-S,(w)ep  +  +  (i  -  l)i?Tpd) 

0 


(197) 


5(w)  -  RT  Pd  +mgRTEz  j  f  a  (l  -  £)  RT  Pd  +RT  Pd 

d>zd 


The  Y1(v)9l  terms  in  (197)  are  now  grouped  to  yield 


d 


ui  =  jt(Bpl)  (~krr  + 


-  “d  -  l)RT^  ~  G(R)  +  RT n -ie, 

Uzd 


~(B  l)kr  (  ~krf  + 


<b: 


-S(u)f.  1 

0 


P  meP 


+  Bu 


z 

0 


-u 


^  y—S(u)ev  +  G(R)  —  RT  Pd  +B\U\ 
0 

H^-S^v  +  ^B^  +  GiR))) 

0 

i  i~B{uj)ep  +  ±ev  +  (£  -  1  )RT'Pd) 
0 


(198) 
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(5-1)  (mRTPd  ~~  aRTPd  -  RT  Pd  +mgRTEz^  +  a  (l  -  B)  RT  Pd  +RT  Pd 

Wzd 

\  Yl(i)0'  1  +  (B-1)  \  i  _  ,  r  sn(«)9 1 1 

,//  '•  '  [  0  ]'<•'[  0  J  '  <*  ;  [  0 


Hence,  the  equation  (198)  can  be  made  more  manageable  in  the  form 


U1=  jt{B^)(t>a^R^^)  +  (b/)Mp,r,v,u)  +  (yjt{B^)(j)c+(B^)(j)d{v)^  Y1{v)9l  (199) 


where  by  the  variables  (pa(p ,  R,  it,  cu),  (j)b(p,  R,  v,  Cj )  G  $i4  have  the  definitions 

(j)a(p,  R,v,lu)  =  -krr  +  U  1, 

4>b(p,  R,  v,  to)  =  —krU2  —  U3  —  U4:  —  U5  —  U6 


(200) 

(201) 


where  U 1  = 


....  \  ~^V  ~  «(ih  -  l)RTPd  ~  G(R)  +  rT  Pd  ~Re% 


U2  =  (  —Ayr 


U  3  = 


-S(u>)rp-  Bep 
0 


U  4  = 


U  5  = 


176  = 


m  — ■S' (d>)e^,  +  B\U\  +  G(R)  —  Rt  Pdj  (202) 

0 

'  H^-S^v  +  BiB^  +  GiR)))  1 

0  J  ’ 

’  hi  (-S,(w)ep  +  Bev  f  (£  -  l)i?Tprf)  ' 

0  J  ’ 

’  S(u)  (tRTPd  ~  «^TIbi  -  RT  Pd  +mgRTEz )  +  a  (l  -  £)  RT  Pd  +RT  Pd  j 


in  which  kr  G  $i4x4  is 

kr  i  0  0  0 

,  _  0  jfcri  0  0 

kr=  0  0  kn  0 

0  0  0  kr2 

and0c(-)  G  $i4x3  and  <f>d(v)  G  $i4x3  are  two  regression  matrices  given  by 


(203) 


0C  =  n  3  ,  (204) 

t^lx3 

4>d  =  (-~h)  J3  I  +  (  kH‘  1  ■  (205) 

V  m  /  |_  t2iX3  J  I  12ix3 

Note  that  the  following  manipulation  was  used  to  facilitate  the  form  of  (j)c  G  Ji4x3  and  <f>d(v)  G  di4x3: 


Y1(v)61 

0 


Y1{v)91  =  BTzYl{v)91  G  d£4. 


Therefore,  the  definition  of  Y3(p,  R,v,lj)93  in  (76)  can  now  be  implemented  using  (199)  and  (193) 
to  produce  the  parameterization  Y:i(p.  R,  v.  uj)  G  J£3x27  and  03  G  T27.  The  final  form  is  given  by 
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Y3(p,  R,  V,  u)  =  [Y31,  Y32,  F33,  F34,  F35,  y36] ,  with  elements 


Y31  = 


Y32  = 


Y33  = 


YM  = 


Yu  = 


Yu  = 


4>b2 

S3 


(jJ1lo3 


—U2t03 
LUlLU3 

£a!  (S±\  _  f  cB_ ' 


^1^2  U)\U)2  (  (^)  _  (§)  0a4  +  ~  %4>b4 


uji  uj\  |a;i|  0 

0  0  CU2 

0  0  0 


S3 


0 

0 

0  0 

|ui| 

0 

0 

^ v2  0 

0 

0  0 

(Pd23 

0 
0 


0  0  0 

0>2  |a>2|  0  0 

0  u)3  a>3  |a^3 1 

0 


0dll 


v\  0 


( _ L  (si.\  1  _  £mi 00  | 

\  S3  \c4>S  '  S3  c<p  J 


Vi 


(206) 


‘Pdll 
'  S3 


V\  |i>i|  0 


<t>d!2 

S3 


V2 


0 

0 


S3 

4>d22 

S3 

0 


( <Pdl2t4>  4,d42°s\  n 

\~3 - — )V 2  U 


^V3  0 


0 

0 


—  <Pd  1 3  J1 

63  3 

0  u3  0 


<t>d23 

S3 

0 


V2  \V2\ 

0 

0 

^32V2  \V2 1 

0 

0 

(  <t>d,12t<t> 

V  s3 

v3  \v3 1 

0 

0 

3  |v3| 

0 

0 

(  4>dl3trf’ 

1  63 

c0  y 


c4>  I 


v2  \V2 


V3  |^3 


and  by  03  —  [031,032,033,034 


CO 

CO 

CO 

with  elements  given  by 

031  — 

Oh 

to 

CO 

01  02 

03 

04  05 

06  ]  , 

032  = 

[  Judi  <72201 

CO 

01102 

02202 

03302 

]. 

CO 

CO 

[  Jlld3  <02203 

03303 

01104 

02204 

-t3 

CO 

]. 

CO 

[  <01105  <02205 

03305 

01106 

02206 

03306 

]• 

(207) 


F  Development  of  Centrifugal/Coriolis  terms 

Consider  the  Centrifugal/Coriolis  force  equations  in  (52),  it  can  be  rewritten  as  follows 

N*  (R,  T,  x)x0  -  N*  (R,  T,  x0)x0  =  [V*  (R,  T,  x)  -  N*  (R,  T,  x0)\  xQ. 

First  we  consider  the  first  term,  N*(R,T,x),  on  the  right-hand  side  of  above  equation 
N*(R,  T,  x)  =  Dt(R,  T)  (C(R,  T)D(R,  T )  -  MD(R,  T,  x)) 


(208) 


where 


67(00,  T)D{R,  T )  -  MD(R ,  T,  x) 
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-mS(uj)  03x3 
Osx  3  S(Juj) 


'  RT  03x3 

mis  03x3 

i(RT)  031 3 

_  03x3  T~\Q) 

03x3  J 

o3x3  i(T-im\ 

mS{uj)RT  03x3 

i 

&l». 

b~i 

03x3 

Osxs  S(Ju)T~\Q) 

03x3 

j&cr-'m . 

-mS{uj)RT  -  mji(RT)  03x3 

o3x3  s{Ju)T-\e)  -  j£{T-\e)) 

03x3  O3X3 


03x3  S(JL0)T~\O)-j£t(T-\O)) 
Then,  (208)  becomes 


N*(R,  T,  x)  = 


R  O3X3 

03,3  T~T(Q) 


O3X3  O3X3 

03x3  SiJ^T-^Q)  -  JiiT-^Q)) 

0^x3  O3X3 

.  o3l3  T-T (0)  [S(j«)r-'(e)  -  J|(T-He))le] 


(209) 


In  a  same  fashion,  we  can  obtain 

]\f*  (  D  T  ™  \  _  ^3x3  03x3 

1  ’  ’  [Oax3  T-T(0)  [SfJwJT-'fe)  -  J|(T-1(e))|e.e„] 

Then,  the  difference  between  (209)  and  (208)  yields 


(210) 


N"(R,T,x )  -  N'(R,T,x„ ) 

O3X3  O3X3 

..  o3„3  r-ye)  [s(jw)r-'(9)  -  J|(r-'(0))]  - r^ie)  [s(jw„)T-'(e)  -  ^(T-'(e))ie,eJ 

03x3  O3X3 

L  oM  r-T(0)  [S(jw)r-‘(e)  -  s(ju,0)T-'m  +  T-ye)  [- J|(r-'(e))  +  J|(r-'(e))|^J 

(211) 

The  first  term  in  the  last  matrix  in  (211)  yields 

r~T(0)  [S'( Jo;)T_1(©)  -  S(Jcv0)T~\Q)] 

=  T~t(0)S(Jlu  -  Jcu0)T-\e) 

=  T~T(6)R(JT-1(6)(6-eo))T-1(0) 

=  T~T(0)S(jr-1(0)s_0)T-1(0)  (212) 

where  from  the  definition  in  (41)  we  can  represent 


n  -  rr»  _  - 

(J  -  O'  o/n  - 


by  defining  as 

sjp  =  p  —  po  G  3?3  and  <s_0  =  0  —  0O  e  Ji3. 
The  second  term  in  the  last  matrix  in  (211)  yields  using  the  (138) 

r~T(e)  +  J^r-ne))!^ 


P~Po 

_A 

S-P 

0  -  0O 

s_0 
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(213) 


where 


=  T-T(0)  '(©))©  + 

=  -T-T(e)j-^(T-\e))  (e  -  e, 
=  -r-T(e)j^(r-1(e))s_e 


|(r-1(e»  =  ^r'(e»eer 

7^7-(T-1(©))  e  S3*3*3  as  a  tensor. 

<70 


Thus,  the  matrix  (211)  yields 

IV*  (i?,  T,x)-N*(R,T,x0) 
03x3 


o. 


3x3 


d 


o3x 3  T-T(e)s(jT-1(e)s.e)T-1(e)-T~T(e)j—(T-\e))s.e 

Hence,  the  Centrifugal/ Coriolis  equations  in  (52)  can  be  defined  as 

N*(R,  T,  x)xa  —  N*(R,  T,  x0)x0 
=  (N*(R,T,x)-N*(R,T,x0))x0 

03x3  0.3x3 

o3x 3  T-T(e)s(jT-\e)s.e)T-\e)  -  t~t(q)j-^(t-1(q))s.q 

Osx3  @3x3 

o3X3  T-T(e)s(jT-\e)s.e)T-1(e)  -  t-t(q)T(j,  0,  *_©) 


Xr, 


o 

•S5, 

o 

o 


3x1 


r-T(0)5(jr^1(0)s_0)r-1(0)0o-T-T(0)r(j,0,s_0)0o 


o 


3x1 


-T-r(0)S(T-1(0)0o).7T-1(0)s_0  -  T-t(Q)Tc(J,  0,  0o)s_0 

03x3  0.3x3 

03x3  T~T(0)  ^-S(T~1(Q)Q0)JT~1(Q)  -  Tc(  J,  0,  0O) 

N*(R,T,x0)s 


S-P 

s_0 


where 


N*(R>  T,  x0)  = 


O 


3x3 


o 


3x3 


03x3  r-r(0)  5,(T-1(0)0O)  JT_1(0)  -  TC(J,  0, 0O) 

and  the  followings  are  utilized 

f(j,e,s_e)  =  j-^(T-\e))s.e, 
zu.  i-i.  =  r(j,e,s_e)e„. 


(214) 

(215) 


(216) 


(217) 

(218) 


G  Figures 
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Tracking  Demonstration 


Figure  3:  Output  Feedback  Tracking  Demonstration 


Position  Tracking  about  y  in  Inertial  Frame 


Figure  4:  Position  and  Yaw  Tracking  at  Each  Axis 
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Position  x  Error  in  Body-fixed  Frame 


Figure  5:  Position  and  Yaw  Tracking  Errors 


Translational  Force  Input  ul  about  lifting 


Figure  6:  Control  Inputs  :  Force  and  Torques 
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Velocity  Observer  Output  about  x  expressed  in  Inertial  Frame 


2  3  4  5  6  7  8  9  10 

Velocity  Observer  Output  about  z  expressed  in  Inertial  Frame 


2  3  4  5  6  7  8  9  10 

Angular  Velocity  Observer  Output  about  roll  expressed  in  Inertial  Frame 


2  3  4  5  6  7  8  9  10 

Angular  Velocity  Observer  Output  about  pitch  expressed  in  Inertial  Frame 


Figure  7:  Velocity  Observer  Output 
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